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ABSTRACT 


In this work the physics of advanced microwave anti-lcIng systems, 
which pre-heat Impinging supercooled water droplets prior to impact, is 
studied by means of a computer simulation and is found to be feasible. 

In order to create a physically realistic simulation, theoretical and 
experimental work was necessary and the results are presented in this 
thesis. 

The behavior of the absorption cross-section for melting ice particles 
is measured by a resonant cavity technique and is found to agree with 
theoretical predictions. Values of the dielectric parameters of 
supercooled water are measured by a similar technique at X ■ 2.82 cm down 
to -I7^C. The hydrodynamic behavior of accelerated water droplets Is 
studied photographically In a wind tunnel. Droplets are found to initially 
deform as oblate spheroids and to eventually become unstable and break up 
in Bessel function modes for large values of acceleration or droplet size. 
This confirms the theory as to the maximum stable droplet size in the 
atmosphere. A computer code which predicts droplet trajectories in an 
arbitrary flow field is written and confirmed experimentally. Finally, 
the above results are consolidated into a simulation to study the heating 
by electromagnetic fields of droplets impinging onto an object such as an 
airfoil. Results indicate that there is sufficient time to heat droplets 
prior to impact for typical parameter values and design curves for such a 
system are presented in the study. 
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CHAPTER 1 
INTRODUCTION 


Aircraft Icing has long been recognized as one of the most 

serious meteorological hazards to flight. Ice formation can result in 

loss of aerodynami c efficiency , control and visibility, as well as 

an increase in aircraft weight and the failure of vital communication 

or instrumentation sources. While a variety of techniques have been 

applied to the icing problem, there are still areas where present 

I 2 

systems fall short. ’ 

The work described In this thesis is centered around understanding 
the physics of advanced concepts in aircraft ice protection which employ 
microwave electromagnetic radiation. While understanding the physics 
of advanced concepts was the fundamental thread which united the work, 
it was often necessary, in the course of the work, to cross Into other 
fields in order to obtain a satisfactory understanding of the physics 
involved. Some of the fields in which experimental or theoretical 
work was necessary are Atmospheric Physics, Experimental and Compu- 
tational Fluid Dynamics and the Interaction of Electromagnetic 
Radiation with Hatter. As a result of the research, contributions 
have been made in each of the above fields enroute to the original 
objective of understanding the physics underlying microwave ice 
prevention systems. 

In Section l.l the Icing problem will be defined and current 
techniques used to deal with it v«MI be discussed. In Section 1.2 
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some of the concepts for microwave Ice protection will be described. 

In Section 1.3 en outline of the following chapters and their relation” 
ship to the thesis will be presented. 

I.l The Icing Problem 

Ice forms on aircraft structures when flight Is conducted through 
areas of supercooled cloud or precipitation droplets. Supercooled 
water droplets, which occur commonly In the atmosphere, exist In a 
metastable state. If some structure, such as an aircraft, comes In 
contact with a supercooled drop, then It will begin to heterogeneously 
nucleate and form Ice. The rate of nucleatlon and subsequent freezing 
depends on the structure and on the temperature of the drop. This 
Is the basic mechanism for Ice formation on aircraft. 

For slightly supercooled drops (”5^C to 0®C), the drops freeze 
slowly after Impact and smooth "clear" Ice Is formed. For colder 
temperatures (-20^C to -5^C) , droplets freeze quickly and an opaque, 

I rregular "rime" Ice Is formed. Below ”20^C, supercooled water 
becomes less common In the atmosphere as homogeneous nucleatlon begins 
to occur and at temperatures below *I|0^C essentially all water Is In 
the Ice phase. ^ Frozen hydrometers do not contribute to the Icing 
problem, as those particles which do strike the structure bounce off. 

The values of droplet and liquid water content applicable to 
U.S. continental cumuli form and stratiform clouds are presented In the 
design criterion for Ice protection certification In Part 25 of the 
U.S. Federal Aviation Regulations (FAR’s).^*^ However, the FAR's 

neglect Icing from supercooled precipitation droplets, such as 
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freezing rain. Table 1 summarizes the range of droplet parameters 

applicable for cloud and precipitation icing. 

While a description of the meteorological conditions which lead to 

Icing is fairly easy to provide, accurate forecasting of these 

conditions is much more difficult. This is due to difficulties in 

predicting whether a cloud will be in the supercooled-water or ice 

phase, which can change quickly with time. For example, if some ice 

phase particles are Introduced into a supercooled cloud, through 

heterogeneous nucleatlon or some other means, then the cloud will 

very rapidly transition to the ice phase, due to the fact that the 

6 

saturation vapor pressure over ice is lower than over water. 

In addition. Icing zones have been observed to be fairly localized 

even when accurately predicted. ^ 

The uncertainty in predicting icing conditions causes forecasters 

to be conservative and to forecast icing conditions anytime the 

potential exists. The conservative nature of icing forecasting has 

two detrimental effects. The most Important Is the "cry wolf 

syndrome, where pilots become accustomed to flying in forecast icing 

conditions with no difficulty and ignore Icing forecasts in more severe 

conditions. This is borne out by the fact that in 92t of the icing 

accidents between 1973 and I977» icing conditions were correctly 
2 

forecast. 

The second detrimental effect of conservative forecasting Is 
economic. Most general aviation aircraft are not certified for flight 
Into "known" icing conditions. While the definition of "known" 
icing is somewhat airbiguous, many operators choose not to fly in 
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Table I. Parameters of the Icing Problem 


i 

! 



Cloud Droplets 

Precipitation Droplets 

Temperature 

0°C 20°C 

6°C -5°C 

Liquid Water Content 

0 -*• I gm/m^ 

0 -► 3 gm/m^ 

Hean Drop Diameter 

10 -*■ AO microns 

1 -► 5 mm 
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forecast tctng conditions. Cancellations due to Invalid forecasts 

are a drain on the aircraft operators and the economy In general. 

The reason for the conservative approach to aircraft icing are 

the Inherent dangers. Some of these are: 

Loss of aerodynamic efficiency 
Increased drag 
Decreased lift 
increased stall speed 
Increased weight 
Loss of control move 'lent 
Engine failure 

Aeroelastic flutter resulting from change in structural mass 
distribution 

Loss of visibility through windshield 

Loss of navigation and communication antenna 

Loss of cockpit instrumentation sources (pitot, static, etc.) 

Loss of onboard radar 

The above icing-related phenomena can occur singly or multiply with 
varying degrees of severity, depending on the icing conditions, the 
aircraft and the Ice protection equipment on board. 

The available Ice protection devices come in two basic 
forms, anti-ice (no ice Is allowed to form), and de-ice (some Ice forms 
and is subsequently removed). The current techniques all have 
advantages and drawbacks. Electrothermal and hot-air anti-ice 
devices are very effective, but require large amounts of energy in 
that they evaporate all iiqpinging water at a cost of 600 cal/gm (2.5 
X 10^^ erg/gm). Freezing-point depressants such as glycol or chemical 
pastes are efficient but are subject to erosion on the leading edge, 
and combustion problems. Pneu-tic de-icing boots are efficient, 
but tend to be unreliable and svbjecc to misuse. 
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From the foregoing, there Is clearly a need for advancement, both 
in the ability to accurately forecast icing conditions and in ice 
protection once icing conditions are encountered. in addition, a 
better understanding of the Icing problem on the part of pilots, as 
well as researchers, will help increase the safety and productivity 
of flight. 

1.2 Concepts for Microwave Ice Protection 

The use of microwave electromagnetic energy for ice protection 
has been proposed by a variety of researchers for both ice detection 
and ice prevention roles. Remote detection of localized icing zones 

g 

was proposed by Atlas In \SSk. Using meteorological radar, zones 
of moderate to high liquid water content above the freezing level 
can be observed in real time. With the advances in ground-based and 
airborne radar, this approach holds great promise. In 1976, 

Magenheim proposed using microwaves to detect the local accumulation 
of ice on airfoils and helicopter rotors by measuring the change in 
impedance, as ice accumulates on the dielectric coating of a surface 

Q 

waveguide located on the leading edge of the airfoil. The 
technique is fairly successful but suffers from anomalous measurements 
when llquie water Is mixed with the ice. 

Hagenheim also proposed using microwave energy in a de-icing 
system. In this scheme, microwave heating at the ice dielectric 
interface of a surface waveguide, located on the leading edge of the 
airfoil, is used to break the adhesion bond of ice. This technique 
was demonstrated in 1976. However, it has no real advantage over 
electrothermal de-icing techniques, which operate on the Identical 
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i principle, and has some disadvantages in terms of complexity and 

efficiency. 

In 1980, Hansman and Hollister proposed using microwave heating 

■ ‘ 1112 
I to prevent the formation (anti-ice) on aircraft structures. ' 

The concept Is to preheat the supercooled water droplets to above 

freezing prior to impact by a microwave field ahead of the airfoil and 

f 

thereby prevent ice from forming. The potential advantages seen for 
a microwave anti-icing system are low power consumption, low maintenance 
and aerodynamic cleanliness, as opposed to other anti-ice techniques. 

Low power consumption is anticipated due to the saving of the latent 
heat fusion (80 cal/gm, 3*35 x 10^ erg/gm) by circumventing the 
water-to-Ice-to-water phase transitions, and to the ability of 
selectively heating supercooled water droplets. The selective heating 
of water droplets is a result of the strong absorption character’stic 
of water in the microwave regime, whereas snow, ice and metal surfaces 
are poor absorbers. This Implies that the wing need net necessarily 
be as hot as with other techniques, which minimizes convective and 
evaporative losses. 

The advantage of keeping the airfoil as cool as possible can 
be clearly seen in Figure l.I, where the heat loss from a wet surface 
exposed to a tengent*al flow Is plotted as a function of surface 
teinperature. ^ The ambient temperature Is - 20 ^C and the tangential 
velocity Is 60 m/sec. The heat flux increases dramatically above 
0 ®C due to conductive and evaporative losses. In order to run the 
airfoil as cool as possible and still prevent runback freezing 
problems, the optimally efficient anti-icing system Is mo!.t likely 


1 



Heat loss from surface (kw/m*^) 
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Figure 1.1 Heat loss from a wetted surface. 

(calculated from NACA TN 2799’’ ) / 

Ambient temr . 20^C / 

Surface velocity 60 m/s / 
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a mtcrowave hybrid. in a hybrid system a microwave leading edge 

system would be combined with either a freezing-point depressant or 

an electrothermal system on the aft airfoil sections. Advanced 

14 

freezing-point depressant pastes have been successful at suppressing 
runback refreezing, and can operate at ambient temperatures, but are 
subject to rapid erosion on the airfoil leading edge. Electrothermal 
systems which operate slightly above freezing are aiso efficient, 
but have a problem initially heating the droplets. These techniques 
are, therefore, well-suited to being combined with microwave pre- 
heating In order to approach optimal efficiency in cases where 
anti- icing is required. 

Power requirements for typical general aviation parameters are 
estimated for the microwave system to be on the order of 100 W for 
propeller anti-icing and 1 WJ for wing anti-icing. An additional 
advantage of the microwave system is that, neglecting circuit losses, 
power Is only consumed when liquid water is present, and thereby 
has the capability to serve as its own detector. 

An example of a possible first-order microwave anti-icing scheme 
employing surface waveguides is shown in Figures 1.2a and 1.2b. In 
Figure 1.2a, a cross-section of the airfoil showing the dielectric 
Inset for the surface wave is presented. The electromagnetic wave 
propagates along the leading edge of the airfoil and is bound to the 
surface waveguide. The electromagnetic field strength characteristi- 
cally decays exponentially away from the waveguide. The waves are 
launched at the wing root and propagate to the tip, where the energy 
not absorbed by the water droplets is collected and recycled back 
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Figure 1.2a Cross-section of an airfoil with 
surface waveguide. 



Figure 1.2b Schematic of the microwave anti-icing 
system employing surface waveguides. 
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inside the wing to the root In a "race track" pattern. It should be 
noted that the launching- retrieving process is not perfectly 
efficient. The presence of liquid water can be measured by a power 
loss as observed through a directional coupler and detector on the 
return path of the circuit. 

The microwave anti-ice system has some distinct and unique 
advantages over other anti-icing systems. There are, however, several 
questions to be answered. The primary question is whether there is 
sufficient time to heat impinging droplets at velocities of 
aeronautical relevance (80 to 200 mph). While there are other 
questions such as runback- re freeze, electromagnetic field optimization 
and circuit efficiency, these become moot if there is not sufficient 
time to heat the droplets. For this reason, the major thrust of 
this thesis will be to endeavor to determine, in a systematic manner, 
whetfier there is sufficient time to heat atmospheric water droplets 
to above freezing prior to impact. 

i.3 Thesis Structure 

The primary goal of this thesis, as discussed In Section i.2, 
is to understand the physics of a supercooled water droplet being 
heated by microwave electromagnetic energy, as It approaches an 
airfoil. This problem is studied by means of a computer simulation 
in Chapter 6, which was designed to be as physically complete as 
possible. In order to provide the physical background for the 
simulation, some preliminary work of a more basic nature was necessary. 

In Chapter 2, the theory of absorption by hydrometers is 
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reviewed. In Chapter 3> the absorption and scattering of mlxed>phase 
Ice/water particles Is studied, both theoretically and experimentally. 
In addition, values of the dielectric parameters for supercooled v.ater 
droplets are experimentally measured. In Chapter h, the deformation 
and stability of droplets subject to an external acceleration are 
studied, both theoretically and In the wind tunnel. This work Is 
Intended to be combined with that of Chapter 2, where deviations from 
sphericity are found to have a pronounced effect on the absorption 
properties of droplets. 

In Chapter 5, a computer simulation of water droplets Is 
presented and experimentally verified. Finally, In Chapter b the 
results of the previous chapters are combined into a computer code, 
which follows the trajectories and heating of droplets as they approach 
an airfoil. In Chapter 7, the results and contributions of the 
work are briefly reviewed by way of a conclusion. 
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CHAPTER 2 

THEORY OF ABSORPTION AND SCAHERING BY HYDROMETERS 


A brief summary of basic scattering theory applicable to atmospheric 
hydrometers is presented in this chapter. In Section 2.1 some 
scattering terminology is defined. The known dielectric properties of 
liquid water and ice are presented in Section 2.2. in Section 2.3 and 
Section 2.4 the approaches of Hie, Rayleigh and Gans for scattering from 
spheres and ellipsoids are briefly discussed. Finally, in Section 2.5 

the problem of a collection of scatterers is considered. 


2.1 Scattering Terminology 

In the classical formulation of the scattering problem, a plane 
wave Is incident on the scattering object. Shown in Figure 2.1 are the 
cases of Interest here. The wave is electromagnetic and propagates In 
the di recti on wi th wave vector k^ ■ The wave is 

polarized wi th the electric field E^ in the e^^ direction. 

The scattering object is of arbitrary shape and is ini>edded in a 
background of uniform dielectric material characterized by the complex 
dielectric constant e^. For most atmospheric applications, is 
taken to be unity, the value for free space. The scattering object has 
a complex dielectric constant e which can vary spatially. The 
background and the scattering object are assumed to be dielectric in 
nature. Therefore, the magnetic permeabilities of the background and 
the object are unity. 
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Figure 2.1 The scattering problem showing the Incident and 
scattered polarization and propagation vectors. 
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The energy flux of the Incident wave 1$ characterized by its 
Poynting vector 


which is aligned in the n^ direction. 

When the incident wave impinges on the object, energy can be lost 
either to scattering or to absorption by the object. These losses are 
characterized by various cross'sections. Cross-sections have units 
of area and, when multiplied by the incident Poynting flux, yield the 
power flow into absorption or scattering. Addition of the absorption 
cros«-section and the scattering cross-section equals the 

o S 

total cross-section 

o. « o + a. (2.2) 

t a s 

which measures the total energy flow out of the wave due to the object. 

In many cases, the scattered energy is not isotropic. The 
differential cross-section, used. When multiplied 

by |S| I It yields the power flux with polarization e through a 
differential of solid angle centered about the direction defined by n. 
The scattering cross-section Is then found by integrating over the 
total solid angle and all polarizations* 
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For radar applications the energy scattered back toward the source 
is of special importance. Toward this end, the back-scatter cross- 
section is defined. It is the scattering cross-section of an 

isotropic source with a constant differential cross-section equal to 
that of the object in the backscatter or negative n^ direction. 




da(fi - -Rq) 

25 


dfidg 


( 2 . 4 ) 


It should be noted that in the literature, a , a , 0 are 

o S I 

sometimes replaced by Q^, and 0^ is sometimes simply 0. 

The above convention has been chosen to avoid confusion with the "quality 
factor" II in later sections. 


2.2 The Dielectric Properties of Water and Ice 

The dielectric properties of a non-magnetlc material can be 
characterized by the real and Imaginary parts of Its complex dielectric 
constant. 


e(f,T) - e'(f,T) - ie"(f,T) (2.5) 

The real part is a measure of the polarization of the material 
subject to an applied electric field. The Imaginary part is a loss 
term measuring the energy transfer from the field to the material. The 
dielectric constant is a function of frequency f and temperature T. 
The frequency is related to the free space wavelength X by the speed 
of light c. 
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f • 


c 

r 


( 2 . 6 ) 


Unless othervftse noted, e will be assumed to be equal to unity for 
background materials of air or free space. 

Values of the dielectric constant for water and Ice have been 
extensively measured. Several good tabulations have been published, 
with those of von Hippel and Ryde^^ being especially useful. 

The temperature dependence of e' and e" for water at several 
wavelengths is plotted in Figures 2.2 and 2.3 from von Hippel 's data. 

The dependence of e' and the loss term e" on T Is clearly not 
negligible, particularly at the lower temperatures. The absence of 
data for supercooled water is distressing, because of the strong 
dependence of e on temperature. The lack of data Is primarily due to 
the difficulty In maintaining supercooled water, in the liquid state, 
under controlled dielectric measuring conditions. However, supercooled 
water Is Important because it occurs commonly In nature. 

For some applications It Is convenient to define the complex 
refractive index m as the square root of e: 


m ■ fT ■ n - lie 


(2.7) 


where n Is the refractive Index and k is the absorption coefficient 
of the material. It will also be convenient to define 


( 2 . 8 ) 
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Figure 2.3 Imaginary part of the 
dielectric constant versus temperature 
(from von Hippei ) 
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Values of n, k, |k and the Imaginary part of -K given by Ryde^^ 

are shotvn for water and Ice In Tables 21 and 2.2. 


2.3 Absorption and Scattering by Dielectric Spheres 

The scattering of electromagnetic waves by a dielectric sphere is 
one of the few scattering problems for which there exists a complete 
analytical solution. The general solution was first published by Hie 
in 1908 . Hie's work was preceded by that of Lord Rayleigh. In I 87 I 
Rayleigh solved the scattering problem in the long wavelength limit as 
part of his explanation for the blue color of the sky,^^*^® Both 
solutions are discussed in the following. 


2.3»l Hie Theory 

1 8 

Hie's analysis of scattering by a dielectric sphere is simple 
in concept although somewhat complex In detail. The general approach 
will be outlined here along with the results. The reader Is referred 
to Stratton's text^^or Hie's original work^® for the details. 

Hie's approach was to write the field of the incident plane wave 
as the sum of vector spherical wave functions centered at the sphere. 

The incident field gives rise to oscillating charges In the sphere which 
produce secondary fields Inside and outside of the sphere. These 
external and internal fields are also expanded In spherical wave 
functions. The expansion coefficients for the secondary fields are 
found by Imposing boundary conditions at the origin and the surface of 
the sphere. 
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Table 2.1 Values of the dielectric parameters fo*' water (taken from 
Ryde, 1?) 


X(cm) 



10 

3.21 

1.24 

0.62 

T(*C) 





20 

8.88 

8.14 

6.15 

4.44 

n 

9.02 

7.80 

5.45 

3.94 

0 

8.99 

7.14 

4.75 

3.45 

• 8 


6.^:8 

4.15 

3.10 

20 

0.63 

2.00 

2.86 

2.59 

10 

0.90 

2.44 

2.90 

2.37 

0 

1.47 

2.89 

2.77 

2.04 

- 8 



2.55 

1.77 

20 

0.928 

0.9275 

0.9193 

0.8926 

CM 

o 

0.9313 

0.9282 

0.9152 

0.8726 

1 f 

0 

0.9340 

0.9300 

0.9033 

0.8312 

- 8 



0.8902 

0.7921 

20 

0.00474 

0.01883 

0.0471 

0.0915 

Im(-K) 

0.00688 

0.9247 

0.0613 

0.1142 

0 

0.01102 

0.0333 

0.0807 

0.1441 

- 8 



0.1036 

0.1713 


i 
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fable 2.2 Values of the dielectric parameters for Ice (taken from 
Ryde, 17) 


n 1.78 At all temperatures 

0 24 X 10‘^ 

-10 k 7.9 X 10 

-20 5.5 X 10'^ 


|Kl 


2 


0 

-10 Im(-K) 
-20 


0.197 At all temperatures. This is 

for ice of ur.it density, the 
value to be used when 0 is 
diameter of melted ice 
particle. 

(Marshall and Gunn, 1952) 

9.6 X 10'^ 

3.2 X 10'^ 

2.2 X 10"^ 


* Refractive index of Ice is independent of wavelength in 
the centimeter band. 
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In order to calculate the cross-sections once the secondary fields 
are known, the Poyntin 9 flux is integrated over a spherical surface, 
concentric with the scatterer to yield the energy flow. Integrating 
the total Poynting vector, which includes the incident and secondary 
fields, will result In the energy flow into the sphere and thus a , 

d 

Integrating the scattered Poynting vector with only the secondary fields 
will yield the scattered energy flow cMay from the sphere or o^. The 
problem is simplified by assuming that the sphere of integration is large 
enough that the asymptotic values of the spherical wave functions can be 
used. 

The foregoing analysis results in the following cross-sections 


't 2 tt 


(-Re) Z (2n + l)(a + b ) 
, n n 

n»l 


2 ce 

a - 4- Z 


s 27T (2n+ l)(|aj"+ [bj 


Oa - Ot - 




Z (-1)" (2n + l)(a^ - b ) 
n-| " " 


(2.9) 

(2.10) 
(2.11) 
(2.12) 


where X is the wavelength in the background medium. The expansion 

coefficients of the exterior secondary field a and b can be 

n n 

thought of as the nth magnetic and electric multiple coefficients. 

They are made up of spherical Bessel functions and depend on properties 

of the sphere such as the complex refractive index m and a where 
27TQ 

® and a Is the radius of the sphere. 



The Mie cross-sections have been calculated for ice and water spheres 
by a variety of investigators for different vaiues of a and 
temperature. *7,22-25 


1 



2.3»2 Rayleigh Theory 

For small spheres or long wavelengths 


a ■ 


27ra 

X 


« 1 


(2.13) 


the simplifying assumptions of Rayleigh can be made. Physically, the 
above constraint implies that the incident field is spatially uniform 
over the sphere. in this situation the exterior secondary fields 
from the sphere can be replaced by those of an oscillating dipole 'p 

I . (|^) .3 fi . .3 

m +2 

where a is the radius of the sphere and is the incident electric 
field. The cross-sections and are calculated using these 
fields. The absorption cross-section o can be calculated from the 

O 

ohmic losses of the oscillating charges required to create the dipole. 

The preceding cross-sections can also be calculated using the Hie 

formalism by expanding the coefficients a^ and b^ in terms of ot. 

By neglecting terms of higher order than all multipoles higher than 

the dipole drop out. Both approaches yield the following cross- 
23 

sections 


3 
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a. ■ 0 , + 0 , 

(2.15) 

t a s 

“s’ 3* “ 1'^' 

(2.16) 

0^ - a’ Im(-K) 

(2.17) 


(2.18) 


Mhere 

Values of |K and Im(-K) are tabulated tn Section 2.2. Noting 
that a « ^ the scattering cross-saction and backscatter cross* 
section have the characteristic 1/X dependence known as Rayleigh's 
law. 

The Rayleigh cross*sections are of extreme importance in the study 
of scattering from atmospheric hydrometers. They are valid over a 
wide range of useful parameters and are the basis of comparison for 
those refinements which attempt to in.iude such additional effects as 
non*spherica1 shape or non-homo gencous material. 


2,k Absorption and Scattering by Ellipsoids 

The problem of scattering by an ellipsoid of revolution was treated 
in the long wavelength limit by Cans in 1912. Cans' approach is 
essentially an extension of Rayleigh's work. He resolves the incident 
electric field into components along thres perpendicular axes, one 
being the axis of rotation of the ellipsoid. Cans then assumes that 
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the electric field components excite dipole oscillations along the 
ellipsoid axes. Scattering phenomena can then be calculated by analogy 
to Rayleigh scattering. 

if the ellipsoid coordinate system is ^ » with C the 

axis of revolution, then the dipole moments given by Gans are 


P5-9E5 

(2.20) 

p ■ g'E 

(2.21) 


(2.22) 

where . 

^ _ V(m^ - 1) 

♦ (w^ - l)P 

(2.23) 

_ V(«^ - 1) 

i|ir + (m^ - 1)P' 

(2.2M 

for oblate spheroids. 


P - im - 2P' - ~ [1 - sin"’el 

(2.25) 

and for prolate spheroids. 


P - liir - 2P' - iiir ^ 

(2.26) 

inhere V is the volume of the ellipsoid and e is 

the eccentricity 

defined as 

e - V 1 - (B/A)* 

(2.27) 
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Prolate (2.28) 

Oblate (2.29) 

with A and B being the major and minor radii of the ellipsoid. 

The Cans backscatter cross>sections for water, ice and snow have 

been calculated for randomly and preferentially oriented ellipsoids by 

27 

Atlas, Kerker, and Hirschfeld. Thei r work indicates that for water 
ellipsoids oriented with major axis along E|, increases roughly 

linearly with axes ratio A/B up to vaiues of order 10. This effect 
is geometrical and has little dependence on wavelength. Ice and snow 
are seen to be only weakly dependent on geometrical factors due to their 
generally small dielectric constant. 

The absorption cross-section for oriented ellipsoids will be of 
particular interest in later sections and was not calculated by Atlas, 
et al . Therefore o will be calculated here using the Gans approach. 

Q 

The Incident niectric field is assumed to propagate perpendicular to 
the axis of rotation of the ellipsoid. The problem then becomes 

essentially two-dimensional. in Figure 2. A the orientation angle 6 

* 

is defined as the angle between and o. 

The absorption cross-section o Is calculated by assuming that 

O 

the C and n dipoles have separate cross-sections and which 

are excited by the appropriate components of the incident field. The 

absorption cross-section is then 

2 2 

0^ • <?r cos 0 + 0 ^ sin o 
a ^ n 


k 2 
V - y AB^ 


ii 2 

V - y A^B 


( 2 . 30 ) 
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Oblate Ellipsoid 



Figure 2.4 Orientation of the incident electric field 

with respect to the ellipsOjd coordinate system. 
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The partial cross-sections and o^ are defined by analogy to 
Rayleigh scattering as 


T 


Im(-g) 


(2.31) 


8it' 




(2.32) 


where X is the free space waveiength and g, g' are the Gans factors. 

The ratio of O- and a for water at 0®C to a. of an 
equi volumetric sphere is plotted versus the axis ratio A/B for oblate 
and prolate spheroids in Figures 2.5 and 2.6. Of interest is the strong 
increase in absorption as the eiiipsoid becomes more eccentric. This 
appears to be valid independent of wavelength, subject to the long 
wavelength limit. 


2.5 Absorption and Scattering by a Coliection of Scatterers 

In many cases of significance, the scattering object is actually a 
collection of smaller scatterers. The absorption cross-section o. for 

o 

such a collection is, to first order, just the sum of all the individual 
cross-sections 


Z 0 


• J 


(2.33) 


as long as the scatterers are not in immediate proximity to one another. 
This can be illustrated by considering two Rayleigh spheres of radius a 
separated by a distance r. The scattered electric field from sphere 1 


[ 
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felt by sphere 2 is of order E| where E| Is the incident 

electric field. It Is clear that for Interparticle spacings greater than 
a few radii, the incident field dominates and the particles absorb 
independently. 

Collective effects can be important for the scattering cross- 

sections. In the case of identical scatterers the differential 

28 

scattering cross-section can be written as 

(f) (2.3M 

a1 1 one 

scatterers scatterer 

where jP(q) is the structure factor used in x-ray crystallography. It 
depends on the distribution of scatterers and is defined as 


expdq • (2.35) 

j ^ 

The vector q is the change in wave vector 


and Xj is the position of the Jth scatterer. 

It Is instructive to consider F(^ for the backscatter cross- 
section ^ ■ 21 cq. For a uniform distribution of fixed identical 
scatterers 

Fizt^) - 0 (2.36) 


This is due to the phase cancellation of backscatter signal from 
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one scatterer by another separated by an Integer number of half wave-r 
lengths along the incident direction. If, however, there is some 
relative motion between the scatterers, then the phase cancellation 
averages out and 

F(2k^) - N (2.37) 


where N is the number of scatterers. The backscatter cross-section 
for meteorological scatterers can then be written as 


»b- 


( 2 . 38 ) 


This effect allows the backscatter cross-section to be related to 
such meteorologically relevant parameters as the liquid water content 


and the rainfall rate. 
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CHAPTER 3 

ABSORPTION AND SCATTERING THROUGH THE ICE-WATER PHASE TRANSITION 

The problem of scattering from melting and freezing hydrometers Is 
considered in this chapter. In Section 3.1 the theory of scattering 
from mixed-phase particles is discussed, along with a model for 
scattering by melting atmospheric particles. In Section 3*2, 
experimental techniques used to measure absorption during the phase 
transitions are presented. Finally, Section 3.3 contains the results 
of the absorption cross-section measurements. 

3.1 Theory of Scattering by Melting and Freezing Hydrometers 

The scattering by water-coated ice spheres has been studied 

29 

theoretically by several investigators. In 1951 Aden and Kerker 

extended Hle's theory to the case of two concentric spheres. In 1952, 
30 

Langleben and Gunn calculated cross-sections for water-coated ice 

spheres using the Aden and Kerker results. Additional cross-sections 

31 32 

have been calculated by B.&..en, Herman and Browning. * Experimental 

measurements of the backscatter cross-sections for large melting Ice 

33 

spheres were made by Atlas et a1 in I960. 

Labrum studied the problem of scattering by two confocal 

3k 

ellipsoids in 1952 in an attempt to explain enhanced reflection from 
the "bright band". The "bri^t band" Is observed on meteorological 
radar as a highly-reflective zone located Just below the freezing level. 
The enhanced reflection in the "bright band" can only be partially 
explained by the dynamics of melting precipitation particles. Labrum 
applied the Gans methodology to calculate backscatter cross -sect Ions. 



ko 


He also made some experimental measurements of backscatter from melting 

35 

nonspherlcal ice particles. Because of experimental difficulties, 
however, Lab rum was only able to get general quail to Jve agreement 
with the theoretical backscatter cross-sections. 

In Section 3.1.1, the above theoretical results for spheres and 
ellipsoids will be summarized. In Section 3.1.2, a model of the melting 
and scattering processes which include the theoretical cross-sections 
will be presented. 


3.l.i Water-Coated Ice Spheres and Ellipsoids 

An example of the absorption, scattering, and backscatter cross- 

30 

sections calculated by Langleben and Gunn is shown in Figure 3>l. 

In this case the wavelength Is 3 cm and the equivalent melted diameter 

D of the sphere is 2.4 mm. The cross sections are normalized by 
eq 

their melted values ^^gj^ed plotted against the mass fraction 

of water f 
w 



(3.1) 


where is the mass of water In the sphere and N| Is the mass of 
Ice. 

The behavior of the cross-sections In Figure 3.1 with f^ Is 
typical of those water-coated Ice spheres with radius less than X. 
The scattering cross-sections increase rapidly with f^ to the value 
for water spheres. The absorption cross-section rises very rapidly 
with f^ to more than twice the melted value between f^ ■ .1 and 
f^ * .2 and then gradually reduces to the melted value. 



Figure 
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The Increased absorption of a partial ly-mel ted Ice sphere can be 
understood physically by considering the melted case. The absorption by 
a melted sphere can be thought of as ohmic losses arising from currents 
which create oscillating mu)Mf>o1es to match the boundary conditions at 
the sphere's surface. If an obstruction, in this case an ice sphere, 
is placed in the center of the water sphere, the path length for the 
currents is Increased and the ohmic losses will similarly Increase. If 
the ice sphere is very small there will be little effect, and if the 
sphere is predominantly Ice, then the excited multiples are weaker and 
much of the current will flow in the low loss Ice, resulting in little 
absorption. 

For scattering from melting ellipsoids of revolution, it will be 
assumed that the melting effects of this section and the Cans shape 
effects of Section l.k can be included multlplicatively. Therefore 


o 


g 

®me I ted 


a 

gans • 
^sphere 


sphere 


( 3 . 2 ) 


Here a is the cross*sectIon to be calculated, -=— 2 is the 

^melted 

ratio of partially- to totally-melted cross-sections for a sphere. The 
g 

ratio . - {$ that of the Gans cross-section to an equl volumetric 
^sphere 

sphere and the cross-section of the equi volumetric water 

sphere. 

Backscatter cross-sections calculated as above agree with those of 
Labrum^^ to at least first order. The above simplification is 
considered acceptable considering thav, for most applications, the 
assumption of an ellipsoid of revolution is Itself only an 
approximation of some more complicated shape. 
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3.1.2 Model of Melting of Scatterers In the Atmosphere 

The approximation of atmospheric scatterers as bodies of 
revolution works well for ice, where the low dielectric constant makes 
shape effect unimportant, and for small water drops where surface tension 
dominates. For melting ice crystals, however, a somewhat more complex 
model must be used. In order to model the melting process, snowflakes 
were observed as they fell onto a plate of warm glass. Five distinct 
phases of the melting process were observed, although not every particle 
exhibited all five phases. 

The five phases are shown schematically in Figure 3*2. Phase 1 
is the ice phase before melting. in Phase 2 the crystal begins to melt 
at its extremities, where the heat transfer is the greatest. Small 
water spheres begin to form at the extremities. In Phase 3 enough 
water has melted to coalesce into a water shell. There is, however, 
still sufficient ice structure in this phase to maintai.n a non-spheri cai 
shape. In Phase k the ice has melted to such a point that the 
structural integrity of the ice is gone and the surface tension of the 
water causes the drop to col lapse into a spherical shape. Phase 5 
is that of the totally-melted sphere. 

The absorption and scattering characteristics of a melting ice 
particle are different in each of the five phases. An example of the 
scattering behavior of a melting crystal is shown ir Figure 3>3» where 
0^ and are plotted against the mass fraction f^. In Phase 1 

the particle scatters and absorbs like an equi volumetric ice sphere. 

In Phase 2 the particle can be approximated as a collection of Rayleigh 
scatterers. The absorption cross-section is, therefore, a linearly- 
increasing function of f^. The backscatter cross-section is 
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PHASE t 



PHASE 2 



PHASE 3 


PHASE 1| 



PHASE 5 


Figure 3.2 Schematic representation of the five 
phases of melting. 



1 

f 

5 


I 


i 

I 


I 


1 



ks 


dependent on the form function F (2k^) defined In Section 2.5. It will, 
In general, be a complicated function depending on the orientation of 
the particle and the locations where melting begins. In many cases 
the crystal will have some symnetry, hexagonal or otherwise, and the 
melting should Initiate symmetrically, which simplifies calculation or 
F. The maximum value of F Is N , where N Is the number of water 
spheres. 

In Phase 3, shape effects, as approximated by Gans scattering, are 
combined with the melting effects to calculate the cross'sectlons as In 
Section 3.1.1. In the example shown In Figure 3>3> the cross*sectlons 
are enhanced by shape effects. Implying that the electric field Is 
oriented at least partly along the major axis. It should be noted that 
the shape effects could also degrade the cross-sections to values lower 
than those of a melting sphere, If the electric field were oriented 
along a minor axis. 

In Phase 4, the particle acts like a water-coated Ice sphere for 
o^ and Ojj. Finally, In Phase 5 when the particle Is totally melted, 
the cross-sections are those of a water sphere. 

Freezing water drops In the atmosphere exhibit much simpler 
behavior than the melting case. The drops are generally In the 
metastable supercooled state prior to melting. When crystallization 
occurs. It does so rapidly and relatively homogeneously. The 
scattering behavior of a freezing drop can be approximated by that of 
a sphere or an ellipsoid with a homogeneous dielectric constant 

_ , ^Ice ^I * ^water \ 

”l 


(3.3) 
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which Is the mass weighted average of the dielectric constants 6| 


ce 


B 4 .^- of the constituents, 
water 


3.2 Experimental Methods for Measuring Absorption During Phase 
Transition 

Relatively few experimental measurements have been made to verify 
the theoretical cross-sections of melting and freezing hydrometers. 
Radar observations of phenomena such as the "bright band" yield some 
Insight, but are subject to many uncertainties, due to factors such as 
precipitation dynamics, coalescence, and general uncertainty as to the 
exact nature of the scatterers. 

Controlled experiments which measure the cross-sections of a 
single scatterer are difficult, due to the very small power changes 
Involved. Those measurements which have been made, by Labrum^^ and 
those by Atlas et al,^^ measured backscatter cross-sections under less- 
than-ldeal conditions. Labrum measured for melting hemispheres 
placed in a waveguide. Atlas et al measured 0^ for very large 
melting spheres suspended by a balloon in the near field of a 
meteorological radar. 

A technique for accurately measuring the absorption cross-sections 
of melting meteorological scale drops using perturbation techniques In 
a resonant microwave cavity has been developed. Section 3*2.1 
describes the technique and the particular cavity used In the 
experiments. In Section 3*2.2 the experimental set-up Is discussed. 
Section 3*3 contains the experimental results obtained using this 
technique. 
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3.2 »1 Perturbation Techniques in a Resonant Cavity 

The technique used for measuring the absorption cross-section 
consists, basically, of measuring the change in quality factor Q of 
a high Q resonant cavity where 


energy stored In cavity g, w M 
power loss In cavity P 


(3.4) 


and (0 ■ 2TTf (3.5) 

U > energy stored in the cavity (3.6) 

P ■ power loss in the cavity (3.7) 


The value of Q, can be related to O’ if the drop is non-magnetlc 
and its dimensions are small compared with the scale length of the 
electric field inside the cavity. Under these assumptions the 
scatterer sees an oscillating electric field of strength E^. The 
power absorbed by the drop P^ can be written using o^ and the 
Poynting vector which would result if the electric field were In free 
space . 


P 


d 




(3.8) 


To relate P^ to the change in Q when a drop Is introduced Into 
a cavity, the case of the empty cavity must first be considered. For 
the empty cavity 


1 . 


(I) 


U 

P 


w 


(3.9) 


^9 


where Q is the empty Q and P is the power dissipated in the 

G W 

walls of the cavity. When the drop is placed in the cavity, assuming 
U does not change, the perturbed value becomes 


0) 


P + P . 
w d 


( 3 . 10 ) 


If the droplet Q is defined as 

then 


^d S, Qe 


(3.n) 


( 3 . 12 ) 


Noting Equations 3.8 through 3>I2, the absorption cross-section can he 
wrt tten as 

<*0 " 

2 

If the cavity is such that U is proportional to E^ , then the 
electric field strength will drop out. 

The cavity chosen for drop measurements was a right circular 
cylinder operating tn the TH^Iq node. This mode was chosen because 
of its fairly uniform electric field at the center of the cavity 
oriented along the cavity axis. The TH^jq mode also has vt y small 
surface currents on the end plates near the axis, allowing holes to be 
drilted in the endplates with minimal effect on the cavity fields. 
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The cavity was designed to resonate at 10.66 GHz. Its diameter D 

was 2.153 cm and its length L was 1.229 cm. Figure 3.1| Is the mode 

36 

chart for right circular cylinders, shewing the operating point for 
the cylinder. There Is good mode separation, with the closest mode 
being the TEj^j at \k.k GHz. Since the Is the lowest mode, 

frequencies lower than the TMq^q resonance are cut off. The cavity 
is constructed out of copper and the empty Q is calculated to be 9000. 
The electric and magnetic fields for the TM^jq mode are 

E,(p,t) - (3. Ill) 

6^{p,t) . -IE„J,(^1H) (3.15) 

where p, <(», and z are the standard cylindrical coordinates. With 

the fields known, it is possible to calculate the stored energy U. 
Noting the periodic time dependence and assuming the unperturbed 
dielectric constant to be unity 

^ " TS? [ I Ej(P.O|^ + |B^(p,t)l^ d vol (3.16) 

■'cavi ty 

E * 

“■-r p’ ♦ Ji^ pp 

Jo 

U - O.OOM (O^L) (3.18) 

2 

Note that U has the desired Eq dependence. For the values 
of D and L described above 



(fHD)/2.54xl0’** f(Hz) D(cm) 
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Figure 3.k Mode chart for right circular 
cylinder showing the operating point 
and the closest competing point. 
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0.048 ^ (cm^) 


( 3 . 20 ) 


The absorption cross“sectlon Is simply proportional to the product of 
f and l/Qjj. 


3.2.2 Experimental Set-Up 

The microwave circuit used to measure cavity Q Is shown in 
Figure 3.5. Power was generated by an X-band sweep oscillator which 
varied freguency linearly over a specified range near the resonant 
frequency. Power was transmitted through a variable attenuator and 
by a commercial wavemeter to provide a frequency reference. The 
cavity was connected to the circuit via a short stub to a coaxial tee. 

The coax coupled magnetically Into the cavity fields by a loop oriented 
In the (fi direction located on one of the end plates. The tee was 
connected to the rest of the circuit by flexible coaxial cable. This 
allowed the entire cavity assembly to be placed In a cold t>ox. The 
cable was terminated at both ends by matched attenuators In order to 
reduce reflections. 

After passing the cavlty« power was measured by a crystal detector. 
The crystal had been calibrated previously against a bolometer. The 
output of the crystal went to a digital scope. The scope was 
triggered by the sweeper and had the capability to store traces on floppy 
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disks for later analysis. This feature was used when measuring time- 
dependent phenomena, such as melting effects. 

in Figure 3.6 a typical scope trace of crystal voitage versus time 
or, equivalently, frequency, is shown. Th#* reference marker and the 
cavity resonance are visible. When the frequency is off resonance, 
no power is absorbed by the cavity and the power to the crystal is at 
a maximum. On resonance the cavity is very absorptive and the crystal 
power is at its minimum. The loaded Q of the cavity, which includes 
coupling effects, is 



r 

a 



where f^ is the resonant frequency and 
maximum of the resonance. The cavity Q 
the unloaded Q which is related to 


r is the full width at half 
discussed previously is 
for a circuit such as this by 


a 


0 




( 3 . 22 ) 


where P and P , are the maximum and the minimum powers 
max mi n 

received by the crystal. 

The actual cavity assembly used in the experiment is shown in 
Figure 3*7. It was machined out of 1.5 inch copper rod. The cavity 
diameter D was 2.153 cm and the cavity length L was 1.229 cm. One 
end plate was removable for cavity cleaning and polishing. There were 
holes drilled on axis in each of the 1.5~cm-thick et d plates with a 
2 nnt radius to provide access for drop insertion. The cavity was 
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Figure 3.7 Resonant Cavity Assembly 
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supported horizontally on a vee block In direct contact with a 
thernKNneter to measure cavity temperature. 

Drops of distilled water were supported by surface tension on a 
quartz fiber with a radius less than .1 mm. An exomple Is shown In 
Figure 3.8. Drops were placed on the fiber by a syringe and dimensions 
were measured In situ by a direct-read measuring microscope. The fiber 
was oriented along the z-ax's of the cavity and ran through the holes 
In the end plates. It was held on axis by plexiglass supports on 
either end of the cavity. The entire cavity assembly was portable so 
that It could be moved In and out of a -20^C cold box without disturbing 
the drop or the microwave circuit. Typical cooling and warming curves 
are shown in Figure 

The resonant frequency and Q of the empty cavl ty were measured. 
The resonant frequency was 10.65 GHz and the Q was 8990, which are 
very close to the design values of 10.66 GHz and 9000. The differences 
are attributed to the end plate holes and Imperfect cavity surfaces. 

It should be noted that the empty Q varied somewhat due to oxidation 
on the interior cavl ty surfaces, resulting from exposure to moisture 
and the thermal cycling of the cavity. in order to correct for this, 
the cavity was polished periodically and the empty 0 was measured 
prior to each drop run. 

3.3 Experimental Results 

The experimental results of absorption measurements for water drops 
are presented In this section. The measured cross-sections for room- 
temperature drops are found to agree with the Rayleigh values In 



Figure 3.8 V/ater droplet supported on a quartz fiber. 
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Section 3>3<l> In Section 3>3>2 some measurements of the dielectric 
constant for supercooled water are presented. Finally, Section 3<3«3 
contains absorption data for spherical and non-spherlcal hydrometers 
during melting and freezing. 

3.3«1 Absorption Cross-Sections for Warm Drops 

Absorption cross -sect Ions were measured for distilled water 
drops at a room temperature of 20^C using the techniques described in 
Section 3>2. The resonant frequency of measurement was 10.64 GHz 
* .01 GHz. Drops with diameters varying from 0.5 mm to 2.0 mm were 
measured. The measured absorption cross-section is plotted against 
drop volume V In Figure 3.10. The straight line Is a plot of the 
Rayleigh absorption cross-sections 

°a * Im(-K) (3.23) 

with Im(-K) taken from Ryde^^ at 20°C to be O.OI 883 . The data agrees 
well with Rayleigh theory for these small essential ly-spherl cal drops. 
The scatter in the data is mainly attributed to uncertainty in 
measurement of V. 

3 . 3.2 Dielectric Parameters of Supercooled Water 

With the confidence in the measurements of Rayleigh absorption 
cross-sections gained from the results of Section 3.3.1 > measurements 
were made of as a function of cavity temperature T In order 
to Infer the dielectric parameter Im(-K) from the Rayleigh 
cross-sections In equation 3.23. Measurements of |k| , which can be 
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linearly related to the shift In the resonant frequency f^ by 
perturbation theory, were also attempted. Measurements of |k|^, 
however, were unsuccessful because the very small changes in |Kp 
with temperature were masked by relatively large changes in the 
resonance frequency due to thermal expansion of the cavity. Lack of 
|k|^ measurement is not considered too important as, in applications, 
|k|^ is generally assumed to be constant with temperature. 

Measurements of the supercooled temperature dependence of lm(>K) 
were made by inserting distilled water drops in a warm cavity and 
placing the cavity assembly in the -20^C cold box. The drops were 
assumed to remain in thermal equi 1 ibri urn wi th the cavity. Microwave 
heating of the drop was neglected because of the low power of the sweeper 
(less than 1 mW) and the low duty cycle on resonance (less than 0.001} . 
Liquid drops were observed at temperatures as low as -17°C before 
crystallization. 

Values of 1/Q^j and cavity temperature T versus time for a 
typical cooling run are shown in Figure 3-11- The drop cooled with 
Increasing absorption to -17°C where it nucleated and the absorption 
dropped to the value for ice. At 16 minutes Into the run there was 
a sudden jump in l/Oj which lasted for 12 minutes. Similar jumps 
were observed on subsequent runs, although occurring at slightly 
different times and temperatures. These anomalous ly-high values of 
\ 1/Qjj are not thought to be physically related to changes in Im(-K), 

but rather to some thermal effect in the microwave circuit. In the 
following, therefore, such values will be omitted. 

Figure 3<I2 is a plot of the observed values of lm(’-K} as 
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calculated by equation 3-22, along with the values of Im(-K) 
calculated from the dielectric data of von Hippel^^ and RydeJ^ The 
measured data is a fairly smooth fit to the previous values and extends 
well Into the supercooled range. The real and imaginary parts of the 
dielectric constant could be calculated from Im(-K) by assuming |k|^ 
is constant with its value for O^C. The above technique could be used 
to measure Im(-K) at even lower temperatures if a very clean cold box 
environment was maintained. Values at other frequencies could be 
measured If a suitable cavity was rx>nstructed. 

3.3«3 Absorption Cross»Sections During Phase Transition 

In Figure 3>H of the previous section, a supercooled drop is 
observed to crystallize. The value of 1/Q^ dropped linearly in time 
from a value for liquid water to that of ice in less than 1 minute. 

This is the behavior predicted in Section 3«1 for a freezing supercooled 
drop. The melting behavior of ice particles is expected to be more 
complicated. 

Observations of the melting cross-sections were made by placing 

an Ice particle onto the fiber and into the cavity while the cavity was 

in thermal equi 1 ibrium wi th the cold box at -20®C. The cavity 

assembly was then removed from the cold box and allowed to warm. A 

typical example of the behavior of l/Q^ with time is shown In 

Figure 3*13 for an ice sohere of equivalent melted diameter D of 

eq 

1.15 mm. The absorption behaved qualitatively as was predicted In 
Section 3«1> At 12:30 minutes the drop began to melt. The value 
of l/Oj quicly rose to a maximum at 16 minutes. The absorption then 
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began to decdy and at 20 minutes reached the melted value, which was less 
than half the peak value. 

In order to quantitatively compare the measurements with theory, 

Oa 

l/Q. must be converted to the normalized cross -sect I on ■ ■ ■ and 

^ °melted 

some relation must be assumed between time and the mass fraction 


^w " 


w 




+ M, 


(3.24) 


To relate f^ to time T , assume that the droplet temperature Is 
zero during melting. Also assume that melting begins at time t^ and 
ends at time t^. The heat loss from the drop is 

^ ■ const. T(t) (3.25) 

where T(t) Is the tlme^dependent cavity temperature and the constant 
is dependent only on the cavity geometry. From the warming curve in 
Figure 3.9> it is assumed that the cavity temperature Increases linearly 
with time near zero. Therefore 

^ ■ const, t (3.26) 


The change In water mass can be related to the heat loss by the 
latent heat of fusion L|^. 


n, I da 

TF-r;: a? 


(J.27) 
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since the total mass M -f H, Is constant and notlnq that M., "0 at 

w 1 • w 

tine t„ 


>L(t) 


. fu-. 


o> - *o> 


(3.28) 


f^(t) » const, (t - tp)' 


(3.29) 


Noting that f ■ 1 at time t 


t^ » then 


f (0 

w 


(t - 


(3.30) 


which Is the required relation between mass fraction and time. 

Figures 3.i^t 3.15 and 3.16 show the experimental and theoretical 

a 

normalized cross-sections as functions of f for ice 

<^a melted 

spheres equiv'lent .nelted diameters D of 1.15, 1.6 and 2.0 mm. 

eq 

The general behavior is In good agreement with theory, although the 
theory underestimates the absorption cross-sections by as much as 25% 
at peak value in the i^orst case. The discrepancy is due to the fact 
that the Aden and Kcrker theory assumes dielectric values for water on a 
wavelength of 3 cm and a temperature of 18**C. The dielectric constant 
at these values is less absorptive than the actual experimental condi- 
tions of 2.8 cm and O^C. The increased absorption, along with slight 
deviations from sphericity, could account for the observed difference. 

Itie absorption cross-sections for melting non-spherical ice 
particles were observed, to check for the multiple phase behavior 
predicted in Section 3.1.2. Shaped ice particles were constructed by 
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placing a small water drop on a fiber and immersing the drop in a 
liquid nitrogen bath. After the drop froze, it was removed from the 
bath and an additional drop was added. The procedure was repeated 
until the required shape was achieved. The major disadvantage of this 
method was the inability to create sharp edges due to surface tension 
smoothing of the water drops. 

Figures 3*17 and 3*18 show examples of ice particles which wtre 
roughly oblate and prolate ellipsoids oriented with their axis of 
revolution along the electric field. For the oblate case In 
Figure 3* 17, four of the five melting phases are visible. The 
absorption cross-section started in Phase I like that of an ice sphere. 
Once melting began, the cross-section skipped Phase 2 because of the 
relatively smooth surface and absorbed in Phase 3 like an oblate 
water-coated ice ellipsoid with axes ratio A/B >1.8. At a mass 

fraction of 0 . 2 , the drop began to collapse into Phase of a water- 

coated ice sphere and at a mass fraction of I .0 the drop absorbed like 
a water sphere in Phase S* 

In Figure 3*18 an example of a prolate ellipsoid is shown (note 

the change in scale). In this case there was sufficient surface 

irregularity to show a non-uniform melting Phase 2 up to a mass fraction 

of approximately 0.2 where v.ne very strong absorption, In Phase 3» of 

an oriented prolate ellipsoid with A/B ■ 2.S prevailed. Phase 3 

tasted up to a mass fraction of 0.8 which is artificially high due 

to the structural support of the quartz fiber. Above f ■ 0.6 the 

w 

drop collapsed into the spherical shape of Phases it and 5. 



PHASE I phase 
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In conclusion, the experimental measurements of the absorption 
cross-sections seem to support the theory for water-coated ice spheres 
and ellipsoids and the model for atmospheric melting. The observed 
cross-sections were found to be greater than or equai to the theoretical 
values for water-coated ice spheres. The greatest differences occurred 
at the peak of the melting curve, where for one case the theoretical 
cross-section underestimated the measured value by ZS%> 
difference is thought to be the result of discrepancies between the 
theoretical and experimental values of the dielectric constant. 
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CHAPTER 4 

HYDRODYNAMICS OF ACCELERATED DROPS 

The absorption and scattering cross-sections of water drops were 
shown, in Chapters 2 and 3> to be strong functions of droplet shape. 

In order to accurately predict cross-sections, therefore, it is necessary 
to have some modei of the expected droplet shape. The Cloud Physics 
community has done a great deal of theoretical and experimental work 
on droplet shape and break-up in the quasi -steady-state case of drops 
falling in the earth's gravitational field.'^ In Section 4.1 this 

work is briefly reviewed and the generalization to accelerated drops 
is made. 

In Section 4.2 experimental techniques, employing high-speed 
photography, for measuring droplet shapes and velocities in the wind 
tunnel, are discussed. In Section 4.3 the results of experiments on 
drop deformation and stability are discussed. 

4.1 Hydrodynamic Theory of Water Drops 

For very small water droplets, where surface tension is the 
dominant surface force, the droplets assume a spherical shape. Section 

4.1.1 discusses deformations from spherical shape as a result of 
hydrodynamic acceleration of the droplets by another fluid. Section 

4.1.2 presents the oscillations which result from the restoring nature 
of the surface tension. Finally, in Section 4.1.3 the instabilities 
which arise when the surface tension is no longer sufficient to main- 
tain the drop are presented. 
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4.1.1 Drop Deformation 

The steady-state shape of a water drop being accelerated 
hydrodynami cal 1y by another fluid, in this case air, can, in principle, 
be found by balancing all forces acting on the surface. These 
forces are: 

Surface tension 

Centrifugal force from internal circulation of the water 
inside the drop 

Aerodynamic force from air flowing around the drop 
Hydrostatic pressure gradient within the drop resulting 
from acceleration 

Investigators have attempted to model droplets by including surface 

38-42 

tension with different combinations of the above forces. The 

most complete, although reasonably complicated, method was that of 

38 

Purppacher and Piter in 1970. They included each of the above 

forces, either analytically or semi-empirical 1y , into Fourier expansion 

coefficients in elevation angle detailing the change in radius of the 

drop. The results agreed well for quiescent freely-falling drops 

with equivalent diameter less than 5 mm with the experimental 

39 

results of Pruppacher and Beard, including such second-order effects 
as the dimple observed on the bottom of large drops. For drops larger 
than S mm, the Pruppacher and Piter model underestimates the deforma- 
tion somewhat. 

if it is sufficient to approximate the deformed drop as an oblate 
ellipsoid, which is certainly adequate for calculations of Cans 
absorption and scattering cross-sections, then the simpler approach of 
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Green provides accuracy equivalent to that of Pruppacher and Piter. 
Green's approach, in 1975, was to neglect flow effects altogether and 
include only the surface tension and hydrostatic forces. Using this 
approach, the diameter of an equi volumetric sphere 0^ can be written 
in terms of the major- to-mi nor axis ratio A/B of the ellipsoid and the 
acceleration a as 


! 


i 


I 

t 

i i 


9 eq 


/4q vi/2 [ (A/3)^ - 2(A/B)^^^ 1] ’^^ 

9Pw 


(4.1) 


where g is the acceleration due to gravity, (7 is the surface 
tension (n.b., a is also used to denote cross-sections), is the 

density of water and a/g is the acceleration of the drop in g units. 
Green originally assumed that the drops were falling at terminal 
velocity and therefore a/g was unity. The more general case for any 
quasi -steady-state acceleration is included here. Equation 4.1 is 
plotted in Figure 4.1 for water at 20^C. 


4.1.2 Drop Oscillations 

When drops are perturbed from their equilibrium shape, surface 

tension acts as a restoring force. This restoring force results in 

drop oscillations. The oscillations can be characterized by a discrete 

set of normal modes. These modes were first identified by Rayleigh, 

6,47 

who identified the allowed frequencies for a sphere as 


* _ f2n(n - 1) (n + 2)a i' 

f - [ T"— 3 I 

n p D 
w eq 


, 1/2 


(4.2) 










'I'he fundamental frequency is plotted against In Figure k.2 

for millimeter-size drops. 

In the fundamental mode the drop is flattened along one axis 
initially. A quarter-cycle later, the drop is spherical. After one- 
half of a cycle, the drop is flattened along an axis perpendicular to 
the original. At three-quarters of a cycle, the drop is again 
spherical, and after one full cycle the drop is back to its original 
shape. 

The validity of the steady-state drop deformations described in 
Section 4.1.1 after a change in acceleration is related to the funda- 
mental frequency in equation 4.3. For time scales less than one 
fundamental period l/f^ the steady-state behavior will not accurately 
predict the transient drop deformation. For time scales of several 
periods the steady-state behavior will be the average deformation in 
time. After many oscillations, viscous damping will cause the 
oscillations to decay and the drop defomatlon will be accurately 
determined by the steady-state solution. 

It is instructive to consider the possible effect on the backscatter 
cross-section of the fundamental oscillation. The change in drop shape 
will cause a change in the backscatter cross-section via the Gans 
theory. The cross-section will oscillate at the frequency of the 
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drop oscfllation. An example of a periodic axis ratio A/6 end the 
resulting backscatter cross-section is shown in Figure 4.3. 

The oscillation of o^ is a potential source of doppler noise in 
meteorological doppler radar. The rT|j oscillation results In 
sidebands shifted by ^'f^ , as is shown in Figure 4.4. This is a 
source of confusion in that f 2 is of the same order as the doppler 
shift from particles moving at meteorological velocities. These 
f'^tse sideband signals could cause a broadening in the velocity 
distribution function observed in conditions, such as turbulence, which 
could excite drop oscillations. 

4.1.3 Instability and Drop Break-Up 

In the previous section It was indicated that surface tension is 
the primary restoring force maintaining droplet equilibrium. When other 
forces, such as the hydrostatic pressure resulting from aerodynamic 
acceleration by another fluid, becomes comparable to the surface 
tension, then the possibility of droplet Instability and break-up 
must be considered. Drop break- p is important to models of atmos- 
pheric precipitation as a limit on the maximum size of drops and, as 

48 

Langmuir pointed out, that the small droplets resulting from break-up 
are a major source of condensation nuclei in warm rain. 

Observations of drop break-up have been made by Mathews and Mason 
and others^^ on drops ^ailing at terminal velocity. Drops were 
stable In their ro^ ^hly-oblate shape until they reached a certain 
critical base diameter, where a small concave depress: ion In the base 
deepened until the resulting bag-like shape eventually burst. This 
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process is known as the "bag break-up mechanism" and Is shown schemati- 
cally in Figure 4.5. Measurements of the maximum stable base diameter 
at terminal velocity vary between 6 and 9 mm depending on factors such 
as turbulence. 

Theoretically, there are two basic approaches to the break-up 
problem. The first is essentially dimensional. It defines the Weber 
number as the ratio of the aerodynamic pressure on the drop to 

the surface tension stress. 


N 


we 



(4.4) 


where is the relative wind velocity, is the density of air 

and A is the diameter of the base of the drop. For values larger 

than some critical Weber number the drops are considered to be unstable. 

49 50 

The critical value of N has been measured to be of order 10. ’ 

we 

The second theoretical approach, which is also in good agreement 

with experimental results, is somewhat more elegant in that it considers 

4T 

the physical mechanism for break-up. In 1964 Komabayasi et al. 
recognized that the bottom surface of a large drop falling at terminal 
velocity was in a metastable state with a heavier fluid (water) being 
supported In a gravitational field by pressure over a lighter fluid 
(air). The bottom surface would, therefore, be subject to gravitational 
or Rayleigh-Taylor instabilities when perturbations occurred at the 
bottom surface. 

Komabayasi et al. assumed that the perturbations would occur as 
linear capillary-gravitational waves and that the critical size for 
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drop break-up was when the base diameter A was equal to half the 
capillary-gravitational wavelength. For typical values of the surface 
tension, the above predicts a critical base diameter of 8.S5 mm for 
freely-fall ing drops. 

45 

The Komabayasi model was further refined by Klett in 1970, who 
recognized that, in general, drops flatten with axial symmetry and, 
therefore, instability should occur when the base diameter A matches 
the resonances of circular gravitational waves on the bottom surface. 
Following the notation of Yih^^ in cylindrical coordinates, the equation 
for the perturbation of the bottom surface of the drop S(t,p <^) is 


S(t,p,(|)) ® C J^(kgP) cos(n(|>) exp(at) 


(4.5) 


where C is a constant, ^ is the initial growth rate of the 

instability and J. is the Bessel function of nth order. The 
n 

eigenvalue k^ is the wave number of the surface wave and is deter- 
mined by the following boundary conditions at the drop edge 


J-(k p) « 0 
n g 


(4.6) 


Y - - 0 


P -j 


(4.7) 


The allowed solutions for the perturbations, therefore, are the 
Bessel funu ’on modes for a circular metrtbrane. Values of the zeroes 


9 

“5“ are shown for different modes in Table 4.1. 


The growth rate for each mode can be found by balancing forces 
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Table 4.1 The zeroes of Bessel functions and their derivatives 
(kgA) , from Abramowitz and Stegun 


k A 
2 


Zero # 

Bessel Function if 

J (x) - 0 
n 

~J (x) 
dx n 

1 

0 

2.405 

0.000 

2 

0 

5.520 

3.832 

1 

1 

3.832 

1.841 

2 

I 

7.016 

5.331 

1 

2 

5.136 

3.054 

2 

2 

8.417 

6.706 

1 

3 

6.380 

4.201 

2 

3 

4.761 

8.015 
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at the lower surface of the drop, 
at a rate a, the growth rate can be written as 


Assuming the drop is being accelerated 
51 


a 


a(k A/2) 
9 


(p^ + P^) (A/2) 

W a 


[- 


a(p.„ - P,)(A/2)' 


w 


- (k A/2) ] (i|.8) 

9 


2 

Negative values of a imply oscillatory solutions for the pertur- 
bation while positive values imply instability. The most unstable 
mode is therefore the one whose eigenvalue k^A/2 maximizes equation 4.8. 

Once instability occurs, the perturbations grow quickly and the 
perturbation function S(t,p,(j)) is no longer valid; however, the general 
characteristics of the Bessel function origins should be observable 
even in severely-perturbed drops. The "bag break-up mechanism" is, 
therefore, considered to be generated from a perturbation consisting 
of the mode which corresponds to the first zero of the axisymmetric 
Bessel function Jq. It should be noted that the mode is not 

necessarily the most unstable, as predicted by equation (4.8). This is 
due to internal flow in the drop, which has been shown to occur with 
axial symmetry. This flow tends to suppress modes which are non- 
symmetric. The Jq is, indeed, the most unstable axisymmetric 
mode predicted by equation 4.8 for those conditions where "bag break-up" 
was observed. 

Inspection of equation 4.8 indicates that for larger values of 
acceleration or droplet sizes, higher-order symmetric modes should 
become unstable and that, if internal droplet circulation can be 
suppressed, then non-symmetric modes should also become unstable. 
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Some of these modes have, indeed, been observed and are presented in 
Section ^.3. 

k.i Experimental Technicues for Wind Tunnel Observations of Droplet 
Shape and Velocity 

In order to validate the extension, to arbitrary acceleration, of 
the theories of drop deformation and instability, experiments were 
carried out to measure these phenomena for various accelerations. Drops 
were observed photographically in a wind tunnel. The small drop sizes 
and short time scales required the application of techniques of high- 
speed stroboscopic photography. The experimental setup and photo- 
graphic techniques will be described here. 

Experiments were carried out in the MIT low-turbulence wind tunnel 
located in Building I7A. The tunnel has cross-sectional dimensions 
of I ft. by I ft. in the test section and is capable of velocities from 
0 to 80 mph. An optically-clear plexiglass test section was construed 
to allow photographic measurements. Tunnel velocity was measured by 
a pitot-static probe just upstream of the location where droplets 
were introduced into the flow. 

The photographic setup is shown schematically In Figure 4.6. 

Drops were photographed using a shadowgraph technique to minimize 
exposure time. The strobe had a maximum intensity on-center of l8 x 
10^ lux at one meter and had a full width at half maximum of either 
.8 or 3 ysec, depending on the intensity setting. The strobe output 
was diffused by a buffed mylar diffuser screen which backlit the drops. 
The photographs were taken by a 35 mm. camera with a "macro" lens 
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Figure 4.6 Schematic diagram of the Photographic Set-Up. 
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capable of 1 to 1 magnification. To provide maximum resolution, an 
extremely-fine-grain aerial -photography film was used (trade name: 

Kodak Technical Pan) with a high-contrast developer (0-19). The 
resolution combined with this setup was measured to better than 5 mm 
for high- contrast objects. 

The strobe was triggered by a signal from a firing circuit. For 
measurements of droplet velocity, two strobes could be fired sequentially 
to provide time of flight information. The delay between strobes could 
be set at 12, 24, 64, S3, or 124 ysec, as was appropriate to the 
velocities to be measured. 

The camera was generally operated with the shutter held open (bulb 
setting) to avoid focal -plane shutter effects. The f-stops were 
determined experimentally and a range of f-stops were shot on each 
photographic run. At the beginning of each roll of film, the alignment 
grid on the diffuser screen, used to help determine camera position, was 
photographed. A millimeter scale was also photographed in the focal 
plane prior to each run to provide length calibration. The photo- 
graphic data was anab'zed by microscopic measurement directly from the 
film. 

In the experiments discussed in the following section, a stream of 
water droplets was injected vertically itito the horizontal flow. The 
setup IS shown in Ffgui ./*. wMch schematically views the injection 
apparatus alonci the photographic axis. The droplet stream was 
generated by the nozzle, shown in Figure 4.7b, which consisted of a 
1.27 mm beveled orifice in a flat plate. The bevel was located on 
the water side of the plate to suppress nozzle instabilities which 
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Figure 4.7a Experimental Set-Up to photograph droolet 
deformation. 
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Figure 4.7b Schematic detail of nozzle used in 
deformation experiment. 
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caused the stream to wander out of the focal plane. The stream was 
shielded from the tunnel flow until the midpoint of the tunnel by a 
1 .9 cm plexiglass injection tube. This was done to avoid flow effects 
near the wai i . 

The nozzle v ^s run with tap water at either 15 or 30 PSI and located 
55.5 cm above the bottom of the injection tube. The velocity of the 
drops at the bottom of the tube was measured to be 10.4 m/sec at 
15 PSI and 15.& m/sec at 30 PSI. The mean equivalent droplet diameters 
were 1.4l and .98 mm, respectively. The drops were not particularly 
quiescent at injection due to the lack of sufficient fall distance 
between the nozzle and the injection point for oscillations, generated 
at the nozzle, to be damped. 

The velocity shear below the injection tube was measured by a 
movable pitot-static probe. The results for all velocities were very 
similar to those shown in Figure 4.8 for a free-stream velocity of 
60 mph. The velocity increases with distance below the injection tube 
until, at 1 cm, the flow velocity is at the free-stream value. 

Droplets were observed between I and 3 cm below the injection tube 
in the photograph region shown in Figure 4.7a. In this region, the 
drops have not had time to acquire appreciable horizontal velocity. 
Droplet acceleration is inferred, therefore, by assuming that the 
horizontal velocity is zero and that the velocity of the relative wind 
that the droplet feels, V ,, is the sum of the injection velocity 
and the free stream tunnel velocity U 





(4.9) 
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rel 


2 



+ U 
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The acceleration a can be calculated from the hydrodynamic force 
equation 


V 


T (7)^ C 


2 '2 


I^Vel' 


(4.10) 


where M. is the mass of the drnn, A is the base diameter and C. 

is the drag coefficient. The drag coefficients are extrapolated from 

the measured values for falling drops and are discussed in more 

detail in Chapter 5. The valir of M ij- obtained from the 

w 

photographically-measured dimensions of each drop by assuming that 
the drops are oblate ellipsoids. 

4.3 Experimental Results 

Drops observed in the photographic region, described In Section 
4.2, 1 to 3 cm below the injection tube, were found to exhibit tir :0 
types of behavior. For smail drops or low tunnel velocities, the 
drops flattened into oblate ellipsoids with the axis of revolution 
aligned with the direction of acceleration resulting from the relative 
flow velocity V , . For larger drops or higher velocities, the 
drops also deformed initially into oblate ell ipsoi ds but the ellipsoids 
became unstable and the drops broke up. Drops of the first kind are 
described In Section 4.3.1 and the unstable drops are described in 
Section 4 . 3 .?. 
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^.3«1 Drop Deformation 

The axis ratios A/B of drops deformed Into oblate ellipsoids 

1 /2 

is plotted against the equivalent accelerated diameter (a/g) 
in Figure A. 9. The deformations are much greater than those which 

were predicted by the steady-state theory of Section 4.1. The reason 
for the discrepancy can be seen by considering that the time the drops 
are in the photographic region is approximately 2 msec. By comparison, 
the fundamental frequency of the mean size drop * 1.41 mm) is 

205 Hz which corresponds to a period of 4.9 msec. The drops are, 
therefore, in the photographic region for less than one-half of a cycle. 
The observed behavior of the drops is, therefore, clearly transient 
in nature rather than steady-state. 

Given the above, some iiteresting conclusions can be drawn from 
the data in Figure 4.9. While there is quite a bit of scatter in the 

data due to the transient behavior of the drops, there seems to be a 

clearly-defined limit to the transient deformations. This is 
indicated in Figure 4.9. This limit is the maximum value of the 
oscillating axis ratio A/B of a drop which has suddenly been 

accelerated at a rate a. The maximum axis ratio is just whbt is 

required to calculate the absorption and scattering cross-sections 
for drops which experience sudden acceleration due to velocity shear, 
such as In turbulence or near an airfoil. 


4.3.2 Instability and Drop Break-Up 

At tunnel velocities above 40 mph, droplets were observed to 
become unstable in the photographic region. due to the high velocity 





Figure k,9 Axis ratio versus equivalent accelerated diameter. 
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shear and the resulting acceleration felt by the droplets, higher~order 

31 

Bessel function perturbations, predicted by Klett were observed along 

with the "bag break-up" mode observed by others. The theoretical 

inodes are discussed in Section 4.1.3. 

At least six different perturbation modes were positively 

Identified. These Bessel function modes are denoted by letters in 

Table 4.2 along with the boundary conditions and the eigenvalues 

k A/2 at the boundary. Photographs of each of the six modes are shown 
9 

in Figure 4.10 where the relative wind is vertical in the photographs. 

The axisymmetric modes A and C correspond to the first and second 
zeroes of the Bessel function at the droplet edge. The B mode 
corresponds to the first nontrivial zero of the derivative of Jq on 
the edge. in Figure 4.11 a plot of the jQ(kp) surface is provided 
to aid in visualization. 

Mode D corresponds to the first zero of the non-symmetric function 
J|(kp)cos <^. A surface plot of this function is provided in 
Figure 4.12. The periodic dependence with azimuthal angle, cos <4 
is clearly visible in the photograph. Modes E and F correspond to 
the J 2 (^P)cos 2<^ and J^(kp)cos functions. They are also most 
clearly discerned by their periodic dependence on <4 Mode E is 
twice periodic in one azimuthal revolution, while Mode F is three times 
periodic in one revolution. 

The instability threshold for each mode can be determined from 
the growth rate equation 4.8. It is 
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Figure 4.10 Photographs of droplets becoming unstable 
In the six observed Bessel function modes. 
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Figure Ii.l2 Computer generated plot of the J-(kp) cos^ surface. 
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where k^A/2 is the eigenvalue for the particular mode, listed in 
Table 4.2. The acceleration can be found from the force equation 
4.10. This yield, a relationship for the minimum unstable base 
diameter A in terms of the relative wind velocity 


2 (: 




-P ) 


,1A 


[k9A/2]''* 


(A. 12) 


Observed values of A and V , are plotted for the three 
asymmetric modes along with the instability threshold (equation 4.12) 
in Figures 4.13> 4.14, and 4.15. The threshold was calculated 
assuming 0^ • 1 and was equal to the mean value corresponding to 
a 0.98 mm diameter drop. There is good agreement with the theoretical 
threshold as all observed drops exceed the critical base diameter 
for their particular mode. It should be noted that the reason that 
the three axisymmetric modes were chosen for analysis was simply 
the prevalence of data for these modes, due to their ease of 
identification. The non^symmetrlc modes which were observed also 
fulfilled the instability condition (equation 4. II) but less data was 
available due to mode identification ambiguities at certain viewing 
angles. 

In conclusion, the steady-state predictions of drop deformation 
were seen to underestimate the observed values of axis ratio A/B. 

A maximum transient value of A/B was determined and is shown In 
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Table 4.2 

Surface functions, 
of observed inodes 

boundary conditions 

, and eigenvalues 

Mode Surface Function Boundary Condition 

Eigenvalue (kA/2) 

A 

Jo(kP) 

Jo't> - » 

2.405 

B 

Jgdcp) 


3.832 

C 

Jo(kP) 

- 0 

5.520 

D 

J^(kp)cos 

J^(^) - 0 

3.832 

E 

J2(kp)cos 2 

- 0 

5.136 

F 

Jj(kp)cos 3 


4.201 


* This boundary 
the Bessel 


condition corresponds to the second zero of 
function 
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Figure 4.3> The Bessel function perturbations giving rise to drop 

Ii5 

Instability predicted by Klett have been observed experimentally. 
The perturbed shapes and the conditions for instability were found 
to ••'gree well with the predicted values^ 
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CHAPTER 5 

DROPLET TRAJECTORIES 


The assessment of the trajectories of droplets driven by another 
fluid Is an important part of many fields such as aircraft Icing, 
cloud physics and combustion physics. A computer code which simulates 
two-dimensional droplet trajectories, either in a freely- flowing fluid 
or near an object such as an airfoil, has been written and is described 
in this chapter. In Section 5.1 the computer code is described. in 
Section 5.2 wind tunnel experiments designed to validate the simulation 
are discussed. In Section 5.3 some results of calculations for trajec- 
tories near an airfoil are presented. The simulation Is extended 
further In Chapter 6, where the microwave heating of droplets near an 
airfoil is considered. 

5.1 Computer Simulation of Droplet Trajectories 

The first attempts to calculate droplet trajectories In flowflelds 
perturbed by objects were made in 19^0. Glauert^^ studied trajec- 
tories near circular cylinders which became the basis of Instrumentation 

5<» 

to measure the liquid water content of clouds. Kantrowltz studied 

trajectories and their relationship to aircraft icing. Both Glauert 

and Kantrowltz integrated the droplet differential equation of motion 

analytically. The approximations required to do this limited the use- 

55 

fulness of the solutions. In 19^5* Langmuir and Blodgett developed 
a more general approach for the trajectories near simple geometriea) 
shapes which required the use of a differential analyzer. 
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The technique of iteratively solving the droplet equations of 

motion, by finite difference techniques, applicable to arbitrary two- 

56 , 

dimensional flow fields, was proposed by Bergrun in iS^?* The 

computations were done by hand, which made calculating even one 

trajectory a major effort. In 1953 e c’-oup led by Brun and Rinaido 

developed an analog computer which mechanically integrated droplet 
57-59 

trajectories. They computed trajectories near cylinders and 

NACA 65,-200 series airfoils. These analog computations, vyhile 
impressive, were limited due to the difficulty in modifying the programmed 
flow fields. 

The calculation of droplet trajectories, by the iterative solution 
of the droplet equation of motion, is weii-suited to the capabilities 
of modern digital computers. Very little work has, however, been done 
in this area. Those codes which do use digital computers have generally 
been developed by aircraft companies privately and are proprietary. 

A computer code to calculate droplet trajectories in an arbitrary 
two-dimensional flow field has been written and is described in this 
section. Section 5.1. 1 discusses the droplet equation of motion. In 
Section 5.1.2 the iteration algorithm is presented with a flow chart of 
the code. Finally, in Section 5.1.3 a discussion of the appropriate 
drag coefficients, to be used for water droplets, is presented. 

5.1.1 Droplet Equations of Motion 

The equation of motion for a droplet being driven by another fluid 
can be written from the hydrodynamic drag equation as 

“ T ^a (o’ - ^ 


(5.1) 
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where a and v are the acceleration and velocity vectors of the drop» 

(f Is the velocity vector of the fluid, is the. mass of the drop, 

p is the density of air, A 4 . is the area of the drop perpendicular 
d 

to the relative wind and Is the drag coefficient of the drop. It 

is noted that 


u - V 


(5.2) 


is the magnitude of the relative wind felt by the drop, as was discussed 
in Chapter h. The drop mass is 




(5.3) 


where Is the diameter of an equi volume trie sphere and p is the 

eq w 

density of water. 

In order to simplify the use of equation 5.1, it is assumed that 
the area of the droplet is the value for an equi volumetric sphere 
and that any deviation from sphericity will be included in the drag 
coefficient C^. Therefore 




(5.W 


Combining equations 5*2, 5 . 3 , and 5.4 with equation 5.1 yields 

• 5 *^d r D ^ 

^w eq 


(5.5) 
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The Reynolds nutnber Re for droplets Is defined as: 


Re - rf 0 


M *^eq re1 


(5.6) 


where y is tho dynamic viscosity of air plotted against temperature 
in Figure 5.1 from the values in the Smithsonian Meteorological Tables.^^ 
The equation of mrtion (equation 5.5) can be written in the simplified 
form as 


- IT ^ 


(5.7) 


where 


K 


a 


1 

Iff 


P 

w 2 
— 0 
y eq 


(5.8) 


which is constant over a particular trajectory. 

Equation 5.7 is the equation of motion used in the trajectory code. 

It has a form similar to the equation of motion first proposed by 

55 ^d^* 

Langmuir and Blodgett. The coefficient ( ) is presumed to be 

a function of the Reynolds number and is discussed in section 5.1.3; 

K. Is constant (assuming 0^ remains constant); and («T - ^ is the 

• e(| 

relative wind felt by the droplets. 
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figure 5.1 Dynamic viscosity versus 
temperature. 

(from Smithsonian Met. Tables 
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5»1»2 I teratlon Algorl thm 

A flow diagram for the iterative computation of two*- dimensional 
droplet trajectories is shown In Figure 5.2. Three groups of input 
parameters are required. They are the flow field, the Initial droplet 
conditions and the atmospheric conditions. The code calculates 
trajectories for an arbitrary known two-dimensional flow field. The 
velocity field u(x), where x is position, must be input either through 
flow equations or by means of tabulated values. If trajectories near 
bodies are being calculated, then the position of the body surfaces must 
be included. 

The initial conditions for the droplets are the equivalent diameter 
*^eq* initial position x^ and the Initial velocity The 

atmospheric temperature and pressure are required to compute air 
density, viscosity, and K . 

d 

Once the initial conditions have been determined, the equation of 
motion (equation 5.7) Is solved at the initial position. The velocity 
V and position x are then incremented one time step At by the first 
order Euler integration equations 

*i+i ■ Vj At + Xj (5.9) 

^i+l " V, (5.10) 

The value of At was taken to be 10 microseconds In most of the 
simulations, although it could be varied as necessary. 

Once the position and velocity have been Incremented, the code 
checks for hydrodynamic stability of the drop (as determined by the 
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stability conditions in Chapter k) and for col lision wi th the body. If 
either of the above conditions occur, then the code kicks out of the 
iteration loop and outputs trajectories to that point. If the drop 
remains stable and no collisions occur, then the code cycles back and 
calculates acceleration at the new droplet position. The code then 
progresses through the loop. Values of droplet velocity and 
acceleration can also be output along with the trajectories. Additional 
forces such as gravity can be included in the acceleration equation. 


S.1.3 Drag Coefficients 

In previous calculations of droplet trajectories, investigators have 

53 Sk 

used a variety of drag coefficients. Glauert and Kantrowitz assumed 
Stokes (viscous) flow around spherical droplets. Langmuir and Blodgett^^ 
observed that at aircraft velocities and meteorological droplet diameters, 
Stokes' law does not hold. They proposed using drag coefficients 
measured for solid spheres. 

In order to check the validity of the Langmuir and Blodgett drag 

CjRe 

coefficients C^, or equivalently ( 2 ^ ) » The code was used to compute 
the terminal fall velocity of drops in the earth's gravitational field. 
The results are shown In Figure 5.3, along with the we II -accepted 
experimental results of Gunn and Klnzer.^^ The Langmuir and Blodgett 
coefficients predict a higher terminal velocity than is observed 
experimentally. The reason for the discrepancy is that the Langmuir 
and Blodgett coefficients neglect such effects as drop deformation and 
the Increased drag due to turbulent eddies downstream of the drop. These 
effects increase and reduce the terminal velocity. 

in order to improve on the Langmuir and Blodgett coefficients, It 
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was decided to Invert the Gunn and Kinzer experimental results to 

determine a set of empirically modified drag coefficients. The modified 
C.Re 

values of ( - ) are plotted against Reynolds number In Figure 5.4 

along with the Langmuir and Blodgett values. The experimental values 

cut off at Re » 3500, due to drop break-up. It is assumed that the 

C.Re 

values of (- yy) can be extrapolated to higher Reynolds numbers. 
Therefore, from Figure 5.4 

C R6 

- 1 .699 X io"^ (Re)’*^^ Re > 3500 (5.11) 

Using the modified drag coefficients, the terminal velocity was 

calculated by the code and Is plotted in Figure 5.5* The values, as 

expected, agree well with the Gunn and Kinzer data for the atmospheric 

conditions under which the experiments were performed. The values of 

I080 mb pressure and 20^C are typical surface conditions. Figure 5.5 

also Includes terminal velocities calculated for typical icing conditions 

of 750 mb and O^C, where the lower air density causes less drag and a 

resultingly higher terminal velocity. The above modified drag 

coefficients have been incorporated into the trajectory code. For 

Reynolds numbers greater than 3500, equation 5.11 Is used. For 

Reynolds numbers less than 3500, the code Interpolates between tabled 
CjRe 

valueb of ( y j p ) . The values which have been used*are listed In 
Table 5.1 for reference. 
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Figure 5*5 
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Table 5.1 Values of (C^ Re/24) as a function 
of Reynolds number. 


Re 

Cj Re /24 

0.00 

1.00 

0.05 

1.009 

0.1 

1 .013 

0.2 

1.037 

0,4 

1.073 

0.6 

1.103 

0.8 

1.142 

1.0 

1.176 

1 .2 

1.201 

1.4 

1.225 

1.6 

1.248 

1.8 

1 .267 

2.0 

1.285 

2.5 

1.332 

3.0 

1.374 

3.5 

1.412 

4.0 

1.447 

5.0 

1.513 

6.0 

1.572 

8.0 

1.678 

10.0 

1.782 

12.0 

l.VOl 

14.0 

2.009 

16.0 

2.109 

18.0 

2.198 

20.0 

2.291 

25.0 

2.489 

30.0 

2.673 

35.0 

2.851 

40.0 

3.013 

50.0 

3.327 


Re 

Cj Re /24 

68.7 

3.664 

98.9 

4.409 

134.0 

5.170 

175.0 

5.943 

220.0 

6.683 

269.0 

7.520 

372.0 

9.480 

483.0 

11.47 

603 . 0 

13.69 

731.0 

16.08 

866.0 

18.66 

1013.0 

21.27 

1164.0 

24.01 

1313.0 

27.03 

1461.0 

30.32 

1613.0 

33.81 

1764.0 

37.36 

1915.0 

41.49 

2066.0 

45.54 

2211.0 

50.12 

2357.0 

54.83 

2500.0 

59.89 

2636.0 

65.24 

2772.0 

71.03 

2905.0 

76.86 

3033.0 

83.41 

3164.0 

87.78 

3293.0 

96.05 

3423.0 

103.69 

3549.0 

111.05 
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5.2 Wind Tunnel Validation of Computer Trajectories 

The computer-generated droplet trajectories have been checked by 
comparison with wind tunnei measurements of droplet veiocity. In 
Section 5*2.1, velocity measurements of droplets injected into a freely- 
flowing wind tunnel are compared with predicted vaiues. In Section 

5.2.2 droplets are observed just upstream of a cylinder where the 
flow field is spatial iy varying but well known. The measured 
velocities are compared with the predicted values and good agreement is 
obta! ned. 

5.2. i Droplets injected into a Uniform Fiow 

In order to verify the computer-predicted trajectories, droplets 
were observed in the set-up shown in Figure 5«6. Droplets were 
injected into a freely- flowing tunnel perpendicular to the flow. 

Droplet velocity components along the tunnel axis were measured 
125 cm downstream of the injection site. The measurements were made 
by the double strobe photographic technique discussed in Section 4.2 
and the wind tunnel was the HIT I ft. x I ft. facility described in 
that section. 

The observed droplet velocities for three values of the tunnel 
free-stream velocity (45 mph, 60 irpi., 75 mph) are shown in Figure 5.7 
along with the computer-calculated values. There Is good agreement 
over the observed 0.15 to 1.0 mm range of equivalent droplet diameters. 
At the observation site the small droplets, due to their low inertia, 
are close to the free-stream velocity, while the larger droplets are 
slower and still accelerating. The discrepancy between observed and 
predicted values of velocity are attributed to difficulty In measuring 
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Figure 5«6 Experimental set up for droplet velocity 
measurements in a uniform flow. 
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O^q, slight variations in the tunnel velocity, and perturbed Initial 
conditions due to droplet break-up or other injection effects. 

S.2.2 Droplet Trajectories Near a Cylinder 

In order to confirm that the comptuer code accurately predicted 
trajectories, in regions of spatially-varying flow, observations of 
droplet velocity just upstream of a cylinder were made in the wind 
tunnel. A cylinder was chosen because the two-dimensional flow ahead 
of a cylinder Is well known to have a simple analytical form. 

An example of the trajectories generated by the code Is shown In 
Figure 5.8. In this case the cylinder had a 10 cm radius and the 
free-stream velocity was 60 m/sec flowing to the left. Four droplets 
with diameters of 5, 10, 20 and 40 microns were Incident on the cylinder 
from a position 1.25 cm above the stagnation streamline. The effects 
of inertia can be clearly seen as the smaller droplets are turned by 
the flow, while the 40-mlcron drop continues on and impacts the 
cylinder. The effects of Inertia on the impingement trajectories will 
be discussed further in Section 5*3. 

The experimental set-up used to measure velocities near the 
cylinder was very similar to that described in Section 5.2.1 and is 
shown in Figure 5.9. The photographic region was centered above the 
stagnation streamline ahead of the cylinder. Trajectories ahead of 
two cylinders of diameters 11.25 and 2.4 cm were observed. The 
cylinders were painted black to minimize stray reflections in the 
photographs. The free-stream velocity was varied from 45 to 75 mph. 
Droplets were injected perpendicular to the flow 125 cm upstream of the 
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Figure 5 >9 Experimental set*up for droplet 
ve’^city measurements ahead of 
a cylinder. 
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cylinder axis and the velocities were measured by the double strobe 
technique. 

Several hundred droplet velocities were observed and compared to 
the computed values. The measured and computed drop velocities were 
found to agree with a standard deviation of less than of the predicted 
velocity for accurately^photographed droplets. 

As an aside, it was observed that when la^ge droplets impacted 
the cylinder, they splashed in rather constant patterns. An example 
is shown in Figure 5.10a. The result of the impact of a large drop 
was a collection of smaller residual droplets with velocities of order 
one-tenth of the free-stream velocity directed radially away from the 
point of first contact with the cylinder. Splash patterns of this 
sort are thought to be the cause of the double-horned Ice build-up 
sketched in Figure 5.10b observed on airfoils under certain icing 
condi tions. 

In conclusion, the experimental evidence supports the results of 
the computer simulations. The code, can, therefore, be applied to 
other situations of Interest with some confidence In the results. 

$.3 Simulation Results Near the Leading Edge of an Airfoil 

In this section, results of droplet trajectory computations near 
an airfoil are presented. The airfoil leading edge is simulated by 
a two-dimensional half body, an example of which is shown in Figure 
5.11. The half body was chosen because it closely resembles the 
front of an airfoil and the velocity field equations have a particularly 
simple form. If the flow Is assumed to be Inviscid, then the velocity 
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Figure 5>10a Photograph of a 
large drop Impacting the cyl- 
inder surface. 
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Figure 5.11 View of the 2D helf body. 
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potential function ♦ Is simply the sum of uniform flow in the x- 

62 

direction and a source at the origin 


♦(x,y) 


U (x + X. 


In 







(5.12) 


where -x^ is the x coordinate of the leading edge 

Xj - h/ir (5.13) 

where h is the half-thickness of the body. The velocity ^ Is 
related to 4 by 

v(xy) « 7$ (x,y) (5.1^) 


The half body equations can be used to simulate droplet impingement 
onto an airfoil. The only weakness Is the assumption of an 0® angle 
of attack a implicit In the half body flow model. Effects such as 
droplet inertia and droplet collection efficiency, which are weak 
functions of oi, can be adequately investigated. In order to study 
affects which depend strongly on a, the particular flow field of 
interest must be inserted in the code. 

In Section 5.3>1( some two-dimensional Impingement trajectories 
are presented. Some results on droplet collection efficiency, and 
their effect on the impinging mass distribution functioivare discussed, 
in Section 5>3*2, the effect of the airfoil on the kinematics of 
droplets is studied. Finally, in Section 5«3«3 some additional 
simulations are presented which demonstrate the flexibility of the model. 
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5«3»l Two-DI mens tonal Impingement Trajectories 

In Figure 5«12, the two-dimensional trajectories are plotted for 
droplets of various equivalent diameter. The effect of droplet inertia 
is apparent, as it was in the case of the cylinders in Section 5.2. 

Small droplets are swept by the airfoil, while the inertia of the large 
drops causes them to resist change in direction and they impact the 
airfoil. in the example in Figure 5*12, the free-stream velocity was 
60 m/sec to the left, and the body thickness was 20 cm. The atmospheric 
conditions were typical icing values of 750 mb pressure (equivalent 
altitude approximately 10,000 feet), ’-IO°C temperature, and gravitational 
acceleration was neglected. In the following, these conditions wil! 
be assumed unless otherwise noted. 

In Figure 5.I3> the tangential trajectory to the half body is 
drawn for a 20-micron droplet. The initial vertical position of the 
tangential trajectory defines the height of the impingement window 
h|(D^^). The impingement window is the zone in which, for a given size, 
all droplets impact the body. For the symmetrical case of the haif body, 
the impingement zone is centered around the stagnation streamline. 
Therefore the height of the impingement window l^| (H^q) Is just twice 
the separation of the tangential trajectory and the stagnation streamline 
well ahead of the body. For non-symmetrical cases, the height of the 
Impingement window is the initial separation between the upper and 
lower tangential trajectories. 

In Figure 5«H, half the height of the impingement window is plotted 
against 0^^ for a 20-cm thick half body. The free-stream velocities 
studied were 40, 60 and 80 m/$ec. It is interesting to note that 
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droplets with D less than 10 microns are always swept by the half 
eq 

body at these velocities. Above 10 microti.., the height of the impinge- 
ment window increases roughly linearly with diameter up to 100 microns 
where the h.(D ) begins to taper to the asymptotic value of 10 cm. 
This linear behavior was also observed by Brun et a1. in their 
simulations. The slight irregularities in the h| curves are not 
physical. They are a remnant of the discrete nature of the sampling 
algorithm and therefore should be ignored. 

A series of trajectories, with varying initial vertical positions, 
are shown in Figures 5.15 and 5.16 for 20 and kO micron droplets. 

The limit of the impingement window is also shown. The 20-m1cron 
droplets are clearly more influenced by the flow than their 40-micron 
counterparts. 

«in^(0e ) 

The differential collection efficiency of the airfoil ° 

eq 

is the ratio of number of dropiets of a given size actually impacting 
the body to the number of drops of that size initially In the 
volume swept out by the body. It can be written In terms of the height 
of the impingement window as 


dn.(D ) 

c _ea_ 


W. 


eq 


H 


(5.15) 


where H is the thickness of the body. The differential collection 

efficiency Is related to the total collection efflcler.cy n by 

c 

integrating over the droplet size distribution f(0_) 

eq 


dTic 

eq 




dO 


eq 


0 


(5.16) 
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Some typical mid-cloud distribution functions are shown in 
Figure $.M. For cumulus clouds, the functions are fairly sharply 
peaked. The mean droplet diameter is of order 20 microns and tends 
to Increase with the size and age of the cloud. Stratus clouds tend 
to have broader distribution functions with similar values of the mean 
droplet diameter. An approximate form of the cloud distribution 
function is the Khrgian-Hazan distribution^ 

f(0 ) ■ const. D ^ exp (3 D /D } (5.17) 

eq eq 6q 6q 


where 0 is the mean effective droplet diameter. 

dH 

c 

The distribution of impinging mass can be calculated from 

eq 

the size distribution function, the differential collection efficiency 
and the volume of the drop. 


dM. 

«<1 


( 5 . 18 ) 


Combining equations 5>17 and 5*18 and noting that H is roughly 
proportional to for cloud-size droplets, and that the volume V 

is proportional to 0^^^ yields: 



const. 



•’V (3 


(5.19) 


An example of the Khrgian-Hazan distribution function and the impinging 
mss distribution function is shown in Figure 5.18 for 0^^ equal 
to 20 microns. It Is interesting to note that the peak value of the 



Figure 5.1 


1 


141 


OF rOJ'.: 





(microns) 


\k 2 


OR!G!T^-- ^ " ;v.,; 
OF POOR QUrtUi ■( 


«<■ 


* *4 


\\ 



c 

<0 


TJ 

9 

O 

u 

c 

♦ 

^ S 

N O 
VI U 

o S 


3 

U 

•M 


lA 

1/1 


o» 

C 

•«p> 

o> 


^ <M 

i ! 



( 

\ 




T 


T 




o 


i 

i 


! ' 


I 


i 


^ ; f ml crons) 



H3 


ORIGINAL PA'i IS 
OF POOR QUALITY 


mass distribution function occurs at kO microns, which Is twice 

the mean diameter D In the cloud. Indeed, the extreme value 

eq 

condition for equation 5*19 yields 





(5.19) 


For raindrops, the differential collection efficiency is approx! 

mately unity. The accepted size distribution function for rain 

63 ik 

Is the Marshall-Palmer distribution, further refined by Atlas. 


f(0 ) « const, exp {-3.67 0 /0^> (5.20) 

eq eq cq 

This results in an Impinging mass distribution of 

dM , 

- const. 0^^ exp {-3.67 } (5.21) 

with the mean value of D for rain being of order several mi 111- 

eq 

meters. 

5 . 3.2 Kinematics of Droplets on the Stagnation Stream, ne 

The general effect of the half body flow field on the kinematics 
of impinging droplets can be observed in the case of droplets which 
flow on the stagnation streamline. Droplets on the stagnation 
streamline are particularly convenient, in that the kinematics are 
essentially one-dimensional. 

In Figure 5.19i the acceleration of droplets is plotted against 
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FIGURE 5. 19 

ACCELERATION VS 01 STANCE ON STAGNATION STREAMLINE 




\ks 

distance ahead of the half bod/, for equivalent diameters of 5, 10, 20 
and kO microns, for the flow conditions listed In Section 5*3.1» The 
acceleration, which is actually a deceleration in that It acts to slow 
the droplets, Is the result of the slowing of the flow due to the pressure 
gradient ahead of the body. Droplets flowing towards the airfoil begin 
to feel some acceleration 30 to 40 cm ahead of the body. The 
acceleration increases until either the larger drops (20 and kO microns) 
impact the body, or the smaller drops (5 and 10 microns) slow to 
approximately the flow velocity a few cm ahead of the body. The 
maximum value of the acceleration is on the order of several thousand 
g units. 

The high values of acceleration that the droplets are subject to, 

result in two major effects. The first is hydrodynamic, In that the 

effective accelerated diameter (— ) D , disc>'ssed in Chapter A, 

9 

a 1/2 

is increased considerably. The value of (~) ' is plotted against 

distance in Figure 5«20, for the droplets described above. The 

peak values of kO to 5S are observed, implying that factors such as 

drop deformation and possibly instability will become important. 

It should be noted that the droplets in the exan^le in Figure 5.20 

a 1/2 

remain below the instability threshold for (— ) ' 0^^ of approximately 

) cm. 

The second effect of the strong acceleration is to increase the 
dwell time of droplets in the region Just ahead of the body. This 
can be seen in Figure 5.21 where the droplet velocity is plotted. The 
droplets are seen to stow appreciably several centimeters ahead of the 
airfoil. The dwell time and, assuming a uniform external droplet 
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field, the droplet concentration, will therefore increase in this 
region. 

it should be noted that the observations made for droplets on the 
stagnation streamline are by no means limited to that special case. 
Droplets in the general region ahead of the airfoil should exhibit 
similar behavior, although the direction of acceleration will be 
variable. In fact, the zone of maximum acceleration is actually located 
slightly off the stagnation streamline Just ahe. J of the body. 

5»3»3 Additional Simulations 

In this section, some additional simulations are presented to 
illustrate the flexibility of the code. The first case of interest 
was that of a large drop impacting the body. As was noted in 
Section 5.2, the large drop splashes into a collection of smaller 
droplets with velocities of order one-tenth of the free-stream velocity 
directed radially away from the point of first impact. 

The splash was simulated by assuming that the point of impact was 
the stagnation point. Simulation droplets were initiated at this 
point with initial velocities of 6 m/sec (one-tenth of the 60 m/sec 
free-stream value). The angle between the Initial velocity and the 
stagnation streamline was varied from 0^ to 90^. Examples of the 
resulting trajectories are shown in Figure 5.22 for the cases of 30°, 

60°, and 90°. For the 60° and 90° trajectories, the droplets are 
clearly swept away by the flow. For the 30° trajectories, the droplets 
return back and Just grazes the body approximately 1 cm along the body. 

The origin of the double-horned icing shape shown in Figure 5«l0b 
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can be Inferred from the simulation results. If the droplets don't 
freeze immediately on impact, most of the residual splash droplets 
will be swept away until some Ice begins to form off the center line. 

The splash droplets, which are moving essentially tangentially to the 
body, will impact at this point and the horn will begin to grow. 

The second horn Is simply a manifestation of the flow symmetry about 
the stagnation streamline. 

Another case of interest Is the effect on droplet trajectories 
of blowing air out the leading edge of the airfoil. Figure 5.23 
shows a series of A0*micron trajectories with and without blowing. 

The blowing velocity at the leading edge is 35 m/sec In this case, 
and the free^stream velocity is 60 m/sec. From Figure 5.23 It is 
clear that It is indeed possible to prevent droplets of a given size 
from Impinging onto an airfoil by blowing. The feasibility of Ice 
prevention by such a technique requires an assessment of the aerodynamic 
penalty, if any, resulting from such blowing and an assessment of the 
power required. 
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CHAPTER 6 

COMPUTER SIMULATIONS OF DROPLET HEATING 

The problem of heating water droplets by microwave radiation in the 
vicinity of an airfoil is studied in this chapter by means of a computer 
simulation. The simuiation is, essentially, an amalgamation of the 
results of the preceding chapters. In Section 6.1 the structure of the 
cof'<e is described. Section 6.2 discusses the model used to maintain 
droplet energy balance in the simulation. Finally, in Section 6.3, 
some results of the droplet heating simulation are presented. 

6.1 Description of the Simulation 

In this section the droplet heating code is described. The code was 
written primarily as a means for obtaining analysis and design criteria 
for aircraft anti-icing schemes which employ microwave preheating of 
water droplets prior to impact. As a consequence, the results of the 
OL.de presented in ths chapter center around droplets which are heated by 
microwave fields as they flow in the vicinity of an airfoil. The code 
is, actually, more general than this specific application and could just 
as easily be applied to problems such as the extinction of solar radiation 
by cloud droplets or the attenuation of radar or communication signals in 
the atmosphere. 

A flow diagram of the code is shown in Figure 6.1. The primary 
structure cf the heating code is the droplet trajectory code discussed in 
Section 5.1. Two sections are added to the trajectory code to account 
for droplet heating. The first accounts for the heating effect of the 
microwave field by calculating the Poynting fiux at the droplet location, 
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given some model of the electromagnetic field and calculating the 
droplet absorption cross-section. 

The cross-section is calculated by first assuming Rayleigh 
absorption (equation 2.17) • The code then checks for drop deformation 
and other non-Ray ieigh effects and corrects the cross-sections by the 
multiplicative factors discussed in Sect^^n 3*1 The following non- 
Rayleigh effects are included in calculating the absorption cross-section 

oa 

large droplets (Hie theory) 
non-spherical droplets (Cans theory) 
phase change effects 

thermal effects on the dielectric constant 
in the above, the magnitude of the drop axis ratio A/B is taken to be 501 
of the maximum transient deformation limit discussed in Section 3<3. 

Given the above, it is possible to calculate the rate of absorption 
of microwave energy by the drop as 



where S is the magnitude of the Poynting vector at the drop location. 

The second additional section in the heating code calculates the 
incremental temperature change over one time step At by considering 
the total heat budget of the drop, including such effects as evaporation, 
melting, diffusive heat transfer and ventilation, along with the 
microwave heating of equation 6.1. The heat budget of a drop Is 
discussed in more detail in Section 6.2. 

Consideration of the droplet heat budget yields the rate of change 
of such droplet parameters as the droplet temperature, T^, the droplet 
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mass, H, and the droplet mass fraction, These parameters are 

Integrated at each time step At by the Euler Integration scheme 
discussed in Section 5.1 


^d,l+l 


^“dt^ * ^d,l 


( 6 . 2 ) 


H 


l+I 


(^)j At + M. 


V,l+1 


df 

(- -) 

' dtM 


At + f 


w,i 


(6.3) 

(6.l») 


The droplet parameters T^, M and are therefore free-running 
parameters In the heating code. 

The initial value of T^ Is generally assumed to be the ambient 

temperature T^. The initial mass fraction and droplet mass are input 

parameters in the code. The droplet mass is represented by the diameter 

of the water sphere of equivalent mass D . 

eq 

1 ■ «-5) 

As In the trajectory code, the droplet parameters and position are 
incremented through the loop until the droplet either Impacts the body, 
misses the body or becomes unstable. 

The code parameters which will be assumed in the following, unless 
otherwise noted, are 

Body ■ half-body 20 cm thick 

Free Stream Velocity ■ 60 m/sec 

Atmosphere - icing conditions (750 mb, -20®C) 
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The electromagnetic field model and the coordinate system used in the 

code are shown in Figure 6.2. The field consisted of a bound surface 

wave wliich propagated along the leading edge of the airfoil perpendicular 

to the flow. The field strength was exponentially decaying with 

2 

distance ahead of the airfoil and had a cos 6 angular dependence 
centered about the stagnation streamline. The field direction is assumed 
to be perpendicular to the radial direction in Figure 6.2. Therefore 

|E| ■ Eq exp{-(r - Xg)/t> cos^ ( j ) (6.6) 

where 


6 • arctan ^ 


(6.7) 


r- lf777 


( 6 . 8 ) 


is a reference angle defining the angular width of the field, 
Eq is the maximum electric field strength at the surface, Z is the 
exponential decay length of the field and x is the x coordinate 

9 

of the leading edge of the half body. 

On the stagnation streamline, the electric field simplifies to 

1^*1 • Eq exp {-(x - Xj)/£) (6.9) 

which is the exponential decay associated with a planar dielectric 
surface waveguide. The exponentiation length Z of such waveguides 
can be controlled by the thickness of the dielectric coating. The 
effect of varying the above electromagnetic field parameters and the 
wavelength X is studied in Section 6 . 3 . 
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Figure 6.2 Coordinate system used for the electromagnetic 
field model. 
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6.2 Droplet Energy Balance 

In order to calculate the rate of change of such droplet parameters 
as temperature, T, droplet mass, K, and mass fraction, f^, It Is necessary 
to observe the heat budget of the drop. The heat budget Is 


dQr 

"It 

•«rf , , -“Wit , -““td 

(6.10) 

where 



«t 

dt 

■ total heat change In drop 

(6.11) 

<*^rf 

dt 

a heat change due to radiation 

(6.12) 

‘‘Qevap 

""it 

■ heat change due to evaporation 

(6.13) 

**^elt 

dt 

• heat change due to melting 

(6.14) 

•«.<. 

TT 

■ heat change due to thermal diffusion 

(6.15) 

The total rate of change of heat Is related to the change In 
droplet temperature by the specific heat at constant pressure 
(Cp » 1 cal/gm^C ■ 4.186 x 10^ erg/gm^C) and the mass N by 



**^t . - „ dT 
"Jt *^p” ^ 

(6.16) 


The rate of change of heat due to radiational heating was written In 
Section 6.1 as 



159 


ORIGIHAL PACS IS 

OF POOR QUALITY 


where o Is the absorption cross-section of the drop and S Is the 

Q 

Poynting flux. 

The change of heat due to droplet evaporation is 


evap , dM . 1 f 

dt 'wv "St ^wv V 'dt - 


( 6 . 18 ) 


where 


■ the specific heat of evaporation of water 
» 600 cal/gm ■ 2.5 x lo'® erg/gm) 


^ ■ the rate of change of mass due to vapor diffusion 


/dHt 

(:rr)n ■ 


the value of jT' for a stationary drop 


the ventilation coefficient which accounts for the enhanced 
vapor diffusion due to air flow around the droplet 


The value of (-^)q may be found by solving the steady-state 
diffusion equation 


3pv -2 
^ ^ ^ Pv 


(6.18) 


where Is the water vapor density and is the vapor diffusion 
coefficient. Calculating the rate of flux of water vapor through some 
closed surface around the drop yields, for spherical droplets^ 


(a?), ■ - ^.d' 


( 6 . 20 ) 


where 0/2 is the radius of the drop and P„ P„ j are the 

w pm VgCl 

ambient and droplet surface vapor densities. If the effect of curvature 
on P^ j can be ignored and the ambient atmospheric is saturater with 


1 ’ 
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water vapor, then P and P . are simply the saturation vapor 

v,a v,Q 

pressures over water and de.'V'nd only on temperature appropriate to the 

pressure to be found. Values of saturation vapor pressure versus 

temperature are plotted In Figure 6.3 from the Smithsonian Meteorological 
60 

Tables. 

Values of the ventilation coefficient 7^, which Is the increase In 
vapor diffusion due to airflow by the drop, have been determined 


empirically by Beard and Pruppacher to be 


f - 0.78 + 0.308 N Re*'^^ 

V sc,v 


( 6 . 21 ) 


where Re Is the Reynolds number and N Is the Schmidt number 

5 Gy V 


sc,v p D 
’ a V 


( 6 . 22 ) 


]i Is the dynamic viscosity of air and p is the density. The value 


of the vapor diffusion Is taken to be 

0 - 0.21 1 


273°K 


(6.23) 


where T and P are the ambient values of temperature and pressure. 

3 3 

The heat loss from the drop can be found by noting that the thermal 
and vapor diffusion equations have Identical form. Therefore, by 
analogy to the vapor diffusion case. 


^ - T, (T, - T,) 


(6.2I|) 


where k is the thermal diffusion coefficient taken to be 

3 

2,k X 10^ erg/cm sec^C in air and 7^ is the heat ventilation 

coefficient, which is taken to be equivalent to 7 following the 

6 

example of Pruppacher. 



Figure 6.3 Saturation vapor pressure over water 


161 




-i - ‘ 


Or TwOH QU*'*-*' 



162 


OR!G:r>J;M i ;3 

OF POOR QUALH / 


The change of heat due to melting of Ice within the drop can be 
change in mass fraction f^ by 
80 cal/gm ■ 3.36 x 10'' erg/ cm) 

(6.25) 


related to the change in mass fraction f^ by the latent heat of fusion 
‘■iw ^'■iw 

*^Snelt _ , .. 


dt 


■-iw” ^ 


In the heating code, no melting is assumed to take place untM the 
droplet is at or above O^C, at which point the temperature remains 
constant and all excess goes into melting ice until f^ becomes unity. 

By considering the above, it is possible to write the rate of 
change equations for the various parameters of interest 

^ - rif '“.S * 7,(L^ 0 ^(p, ^ „) * k,(T, - Tjj) ] 

P (6.26) 


"d“° 

df 


“TT ■ ^^(o /2) f (L D (p - P j) + k (T - T.))] 

dt Lj^ a eq V wv v' v,a v,d a a d ^ 

(6.27) 


dH 


-TT “4Tt(D /2)f L D(p "P ) 
dt eq V wv v'^v,a '^v,a 


( 6 . 28 ) 


in order to check the droplet parameter section of the code, the 
cooling and evaporation of a 1 mm droplet was simulated. The droplet 
was given an initial temperature of 20^C in an ambient saturated 
atmosphere of 0°C. The two cases of a stationary, unventilated drop 
and a freely- fal 1 1 ng, ventilated drop were simulated. The results 
are shown in Figures 6.4 and 6.5. The ventilated drop is seen to 
equilibrate in several seconds, while the unventllated drop requires 
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several tens of seconds. This Is generally consistent with the 

67 

experimental results of Kinzer and Gunn, who found that a 1.35 mm 
drop cooling 7.5^C in subsaturated air took 4.4 seconds to equilibrate. 
There Is, therefore, some confidence that the code reasonably 
approximates the changes in droplet parameters. 

6.3 Simulation Results 

In Figure 6.6, droplet temperature is piuiceJ versus distance 

ahead of a 20>cm'thlck half'body, for droplets impinging along the 

stagnation streamline. In this simulation, the wavelength X was 

1 cm, the maximum field strength was 7.5 kV/cm, the exponentiation 

length of the field ^ was 10 cm, and the free-stream velocity was 

60 m/sec. Droplets of 20, 40, 100 and 1000 micron diameters are 

seen to heat from the ambient temperature of -20^C to values greater than 

O^C upon impact. It should be noted that all the thermodynamic losses 

such as evaporation and diffusion have been Included and that the 

maximum electric field is one-half of the air breakdown value of 

68 

15 kV/cm. It is clear, therefore, that for the above c'>nditions 
droplets can be heated to above freezing before impacting the airfoil. 

Further inspection of Figure 6.6 shows that the smallest drop 
(20 microns) Is warmed the least. This is expected even though the 
smaller droplets dwell somewhat longer In the high field ret ion due to 
their low Inertia. This is due to the small drops being poor absorbers 
by comparison to the larger drops. The trade-off between dwell time 
and absorption cross-section Is what causes the bunching of the 20, 40 
and 100 micron temperatures In Figure 6.6, while for the 1000 micron 
drop the strong non-Raylelgh absorption effects i<imlnate. 
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The effect of varying the electric field parameters is shown for 20 
micron droplets impinging on the stagnation streamline. The 20 micron drop* 
lets have been chosen as a worst-case design point, in that diameters less 
than 20 microns were shown in Section 5.3 to have a low collection 
efficiency and their heating is, therefore, less critical. 

In Figure 6.7 the droplet temperature at impact is plotted versus 
the maximum electric field strength for a wavelength of I cm and several 
values of the exponentiation length. For field strengths above 2.5 kV 
the final droplet temperature increases linearly with electric field 
at a rate of approximately 4.2°C/kV. The exponentiation length is 
seen to have a weak effect on the final droplet temperature, as the 
most effective heating occurs close to the airfoil. Since the ability 
to heat drops does not depend on long exponentiation lengths, then it 
is clearly advantageous to keep the electric field closely bound to the 
airfoil in order to minimize the heating of droplets which miss the 
airfoil. 

In Figure 6.8, the droplet temperature at impact is plotted against 
wavelength for several values of the electric field and a 10-cm 
exponentiation length. The final temperature increases quickly with 
decreasing wavelength for X less than 3 cm. Final temperatures ai r 
freezing are obtainable for X less than 2 cm. 

The field strength required to heat droplets is seen in Figure 6.7 
to have a maximum value on the order of 7.5 kV/cm near the airfoil 
surface with X ■ 1 cm. The power flux for a plane wave of 7.5 kV/cm 
can be calculated from the Poynting vector 


c 

w 


t 


txt 


(6.29) 
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FIGURE 6.7 
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to be 7k kW/cm. This high peak power density implies large power flow in 
the microwave circuit. While this high power is not lost from the system, 
as it recirculates, it creates power handling problems. It may therefore 
be advantageous to keep the fields tightly bound and to operate at shorter 
wavelengths where the maximum electric fields are reduced, as is implied 
by Figure 6.8. It should also be noted that the electromagnetic fields 
are not necessarily plane waves and it may be possible to create a field 
structure which maximizes the electric field while minimizing the magnetic 
field. This would reduce the circuit power flux from equation 6.29, while 
still heating the droplets due to the non-magnet ic nature of water 
droplet absorption. 

In Figures 6.9 and 6.10, the effect of the ambient temperature and 
the flow velocity are plotted for 20 micron and 100 micron droplets. 

The maximum electric field is 7.5 kV/cm with X » 1 cm on a 10 cm 
exponentiation length. In Figure 6.9 with a free-stream velocity of 
60 m/sec, the final droplet temperature decreases with ambient 
temperature, as would be expected. It is noted that the 20 micron 
drop impacts at above freezing, even for ambient temperatures of -30**C 
and that the final temperature of the 100 micron drop averages 
approximately k^ warmer than the 20 micron drop. 

The effect of the free-stream velocity, at a fixed ambient tempera- 
ture of -20®C, is shewn in Figure 6.10. As expected, the high 
velocities provide less time to heat the drops and their final 
temperatures are correspondingly lower. The 4®C to 5®C difference 
between the 20 and 100 micron drops Is again observed. 
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In Figure 6.11 the mass fraction Is plotted versus distance ahead 
of the body. The value of the electric field parameters Is the same 
as above, and the free-stream velocity is 60 m/sec. The Initial value 
of f^ is 0.1. The large (1000 micron) drop Is totally melted by the 
time It is within 6 cm of the body, while the smaller drops are still 
partially frozen upon impact. This is not considered a problem, in 
that most mixed-phase particles in the atmosphere are precipitation 
particles with dimensions on the order of millimeters. 

An example of droplets which impinge from positions oth^r than 
the stagnation streamline is shewn in Figure 6.12. Droplet temperature 
at the nose of the airfoil (x * x^} is plotted against initial separation 
from the stagnation streamline. The droplet size Is 20 mi'rons and 
the reference angle 9 q is 45^. The limit of the impingement window 
is also shown. The plot is somewhat misleading, in that there wit) be 
some additional heating for drops which impact behind the nose of the 
airfoil and the final temperature also depends on diffusive losses. It is 
clear, however, that the droplets which are heated the most are those 
within the impingement limit, and that by judicious choice of the electric 
field the selective heating of only the impacting droplets could be 
maximized. 

In conclusion, a computer simulation has been made of droplets 
being heated by microwave radiation as they impinge onto an airfoil. 

The simulation utilizes the results of the preceding chapters. The 
results indicate that it is indeed feasible to heat droplets to above 
freezing prior to impact under most icing conditions. The effect of 
various electromagnetic field and flow parameters on droplet heating 
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has been calculated, and the results are presented as an aid to the 
design and analysis of anti-icing systems which might employ the 
techniques of microwave heating. The results indicate that wavelengths 
less than 2 cm and maximum field strengths greater than 7 kV/cm are 
required to heat droplets to above freezing prior to impact. While the 

required field strength is high, it is less than the breakdown value of 
air, and the simulation indicates that the fields can be kept quite 
closely bound to the airfoil with little degradation in performance. 
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CHAPTER 7 
CONCLUSIONS 


The primary goal of this thesis was to understand the physics 
involved in advanced microwave anti-icing systems and to determine 
their feasibility in terms of the ability to heat supercooled water 
droplets to above freezing prior to impact. Correspondingly, this 
question was answered by means of a computer simuiation in Chapter 6. 
The conclusion is that it is indeed possibie to pre-heat water 
droplets to above freezing under anticipated icing conditions, as long 
as the maximum electric field exceeds 7 kw/cm and the wavelength is 
less than 2 cm. Detailed design curves for the electric field 
parameters are provided in Chapter 6. 

The work was successful in demonstrating its original goal, 
i.e. that supercooled water droplets can be heated to above freezing 
ahead of the airfoil. In order to create the physically-realistic 
simulation from which the above conclusions were drawn, some 
preliminary work was necessary. The original techniques, methodology 
and results which were developed are most likely as significant as 
the original goal. These results will be reviewed briefly by 
chapter. 

In Chapter 3, a technique was developed to accurately measure 
the absorption cross-section of a single droplet by its effect on 
the Q of a low loss resonant cavity. Using this technique, the 
theory of absorption by water-coated ice spheres of Aden and Kerker 
was experimentally confirmed for the first time. In addition, an 
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extension of the Aden and Kerker theory to melting of irregularly- 
shaped particles was mads and confirmed experimentally. Finally, 
in Chapter 3> an extension of the resonant cavity technique allowed 
the dielectric properties of supercooled water to be measured, at 
X > 2.82 cm, down to -17^C. This new data is shown in Figure 3«12. 

In Chapter 4, the hydrodynamics of accelerated water droplets 
were observed in the wind tunnel by means of high-speed strobe 
photography. The transient deformation limit of the major- to-mi nor 
axis ratio was experimentally determined for droplets suddenly 

accelerated by another fluid. The result is shown in Figure 4.9. 

4c 

In addition, the theoretical prediction of Kiett as to the form of 
Rayleigh-Taylor instability In water droplets was unambiguously 
confirmed as the first six Bessel function instability modes were 
identified photographically. This experimental evidence verifies 
the theory of the maximum stable drop size in the atmosphere. The 
unstable break-up of large water droplets causes the truncation at 
approximately 5 mm observed in the atmospheric exponential droplet 
size distribution of Marshall and Palmer. 

In Chapter 5, a computer code was written to calculate droplet 
trajectories, in any given flow field, and was experimentally 
verified. The code Is useful in analyzing the susceptibility of 
airfoils to icing and to analyzing new ice prevention concepts. With 
the assistance of the code it was determined that, for cloud droplets, 
the larger drops contribute appreciably more to icing due to their 
increased collection efficiency and their increased mass. In 
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addition, the mechanism of the "double-horned** icing shapes resulting 
from large splashing droplets was observed and modeied. 

Finally, in Chapter 6 the above results were amalgamated int^ a 
computer simuiotion, where it was seen that it is indeed possible to 
heat droplets to above freezing, even when evaporative and convective 
losses are included. From the design curves of Chapter 6, some 
conclusions can be drawn as to the practical requirements of a 
microwave anti-icing system. In Figure 6.8 where the strong 
influence of the wavelength on droplet heatirg is displayed, it is clear 
that the minimal wavelength is optimal in terms of heating and that 
a wavelength of less than 2 cm is required. Balancing the desira- 
bility of short wavelengths is the fact that microwave equipment costs 
tend to increase monotoni cal ly with decreasing wavelength. These 
factors indicate that a practical anti-iefng system would operate 
with a wavelength of approximately 1 cm. 

In Chapter 6, for X « 1 cm it was found that a relatively 

high field strength of & roximately 7*5 kV was required to anti-ice 
in all anticipated icing conditions. This field strength is one-half 
the atmospheric breakdown limit. Even though the field strength was 
required to be large. In the results of Chapter 6 it was found t ..t 
it need not be extensive, as a tightly-boun-* it ■ 10 cm) surface 

wave was found to perform as well as a loosely-bound wave it ■ 40 cm). 
This resu' indicates that it should be possible to be quite efficiert 
in only heating the impacting water droplets. With the criteria for 
droplet heating determined by the above results, the practicality 
of microwave anti-icing systems depend solely on engineering 
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efficiency considerations and the ability of the microwave system 
to be integrated with other systems. In a hybrid design, as discussed 
In Chapter 1, to solve the runback- re freeze problem. 
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LIST OF SYMBOLS 
drop major diameter 

drop area perpendicular to the relative wind 

ellipsoid axis ratio 

sphere radius 

drop acceleration vector 

magnetic and electric expansion coefficients of the 
scattered wave 

drop minor diameter 
magnetic field vector 
incident magnetic field vector 
drag coefficient 

specific heat of water at constant pressure 

speed of light 

cylindrical cavity diameter 

equi volumetric sphere diameter 

mean value of 0^ 
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Incident electric field vector 
maximum electric field strength 
electric field strength in mode 

eccentricity of an ellipsoid 
polarization vector 
structure factor 
frequency 

cavi ty resonant frequency 
Doppler frequency 

drop oscillation frequency of the n^^ mode 
mass fraction of water 
average heat ventilation coefficient 
average vapor ventilation coefficient 
Gans factors 

acceleration due to gravity 
airfoil thickness 

magnetic field strength in mode 

body half thickness 
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K 

K 


a 


K 



H 


M. 

I 



impingement window height 
n^^ Bessel function 
dielectric parameter 

constant in the droplet equation of motion 
wave vector 

thermal diffusion coefficient 

gravitational wave number 
i'^cident wave vector 

scattered wave vector 

cavity length 

latent heat of fusion 

latent heat of evaporation 

exponential l&ngch of the electric field 

drop mass 

mass of ice 

mass of water 

impinging mass distribution 


m 


complex refractive index 



I8<i 




we 


sc 


a 

p 

p 

a 

P, P' 

P. 


P P 
min* max 


q 

R 

Re 


number of scatterers 
Webber number 
Schmidt nun^er 
mode numbe r 

real part of the refractive Index 

unit vector in the direction of propagation 

power 

ambient power 

geometrical factors for ellipsoid 

power dissipated by a drop 

power dissipated by the cavity walls 

Hiinimum and maximum power received 

quality factor 

heat 

rate of change of heat 
change in wave vector 
rainfaii rate 
Reynolds number 
radial coordinate 


t f 


I 


i 
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S 


T(t) 

T. 


U 

U 


rel 


X, y, 2 


interparticle spacing 
surface function 
Poynting vector 
incident Poynting vector 
temperature 
ambient temperature 
cavity temperature 
droplet temperature 
time 

time melting begins 
time melting ends 
energy stored in the cavity 
free stream velocity 
fluid velocity 
droplet volume 
relative wind velocity 
droplet velocity 
Cartesian coordinates 


j 
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a 

r 

e 

e' 


X 

X 


y 


c. n, X 


position of the scatterer 
half body scale length 
2ira/X 

instability growth rate 
angle of attack 
full width at half maximum 
complex dielectric constant 
real part of e 

negative of the Imaginary part of e 
airfoil collection efficiency 
angular coordinate 
heating code reference angle 
negative of the imagir.''>‘y part of m 
wavelength 

Marshal I'Palmer coefficient 
dynamic viscosity 
ellipsoid axis of revolution 
ellipsoid coordinates 
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P 




P, 4>» 2 


do 

(2 


dens I ty 

densi ty of Ice 

density of water 

density of water vapor 

cylindrical coordinates 

surface tension 

cross-section 

backscatter cross-section 

absorption cross-section 

scattering cross-section 

total cross-section 
differtial cross-section 

velocity potential function 

solid angle 
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